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are presented to validate the model. 

An analysis of two-dimensional recursive filters is 
presented. A three-dimensional recursive filter is developed 
which exploits the spatial as well as the temporal image 


redundancies. 
A class of hybrid filters is proposed which improves 
the performance of the recursive filters. Several experiments 


with pictures are presented to show the ability of the 
hybrid filters in picture restoration. 

A detector is developed for purposes of target extraction 
from cluttered background images. The detection is independent 
of the target shape. 

A Simulation of the target detection and tracking problem 
1s presented. The target is tracked from frame to frame by 
means of a conventional Kalman filter, which uses the image 
filter as the measurement device. 
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ABSTRACT 


The restoration of images and the enhancement and 
detection of targets in cluttered background are the subjects 
of this research. The statistical approach is used in order 
to exploit temporal as well as spatial image redundancies. 

The images are modeled as a homogeneous random field. 
An autocorrelation function and a method of parameter 
identification are proposed. Experiments with several 
pictures are presented to validate the model. 

An analysis of two-dimensional recursive filters is 
presented. A three-dimensional recursive filter is developed 
which exploits the spatial as well as the temporal image 
redundancies. 

A class of hybrid filters is proposed which improves 
the performance of the recursive filters. Several experi- 
ments with pictures are presented to show the ability of 
the hybrid filters in picture restoration. 

A detector is developed for purposes of target extrac- 
tion from cluttered background images. The detection is 
independent of the target shape. 

A simulation of the target detection and tracking 
problem is presented. The target is tracked from frame to 
frame by means of a conventional Kalman filter, which uses 


the image filter as the measurement device. 
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I. INTRODUCTION 


A. PROBLEM DEFINITION 

This research involves two aspects of digital image 
processing. One is the restoration of images degraded by 
white Gaussian observation noise. The other isthe enhance- 
ment and detection of targets immersed in cluttered back- 
ground images. The statistical approach is used in order 
to exploit temporal as well as spatial image redundancies. 
The emphasis is directed towards the design of recursive 
and hybrid Bayesian filters. 

In the recent past considerable attention has been 
devoted to the application of Kalman filtering to smoothing 
of observation noise in image data. Habibi [4] first sug- 
gested a two-dimensional recursive image estimator as an 
extension of the one-dimensional Kalman filter. Two other 
Similar extensions were proposed [3], [5]. It has been 
shown [10] that these filters do not preserve the optimality 
of the one-dimensional Kalman filter. The most complete 
study of optimal two-dimensional Kalman filtering has been 
performed by Woods and Radewan [11]. They point out that the 
generalization of the Kalman filter to two dimensions can 
be done optimally with an extremely high dimensional state 
vector, which has dimensions on the order of MN (M = order 
of the filter, N = width of the image). Panda and Kak [9] 


succeeded in deriving a vector dynamic model that generates 





the same random field of Habibi's model. Since this model 
is recurSive in only one index, the one-dimensional Kalman 
filter is applied and, of course, optimality is preserved. 
Again, as indicated by Woods and Radewan [11], an extremely 
high dimensional state vector is used, which has the same 
dimension as the width of the image. 

Several techniques have been investigated for background 
clutter suppression and target enhancement. Statistical 
non-recursive spatial filters [6] are suggested for background 
clutter suppression and enhancement of targets of known shapes. 
Two-dimensional Kalman filtering [5] is also proposed for 
the case of targets with arbitrary shapes. These filters 
assume that the image gray level can be decomposed into three 
additive components; target, background and white observa- 
tion noise. The statistical differences of these three 
components are used in order to extract the target, which 


is the desired information. 


B. RESEARCH OBJECTIVES 
In the context of picture restoration and target detection 

as defined above, several specific research objectives have 
heen identified. 

a) Since the model of images plays a fundamental role 
in the statistical approach of image processing, several 
experiments with real life pictures will be accomplished 
in order to validate existing models [1-2] and to suggest 


new ones. 
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b) Generalization of the one-dimensional Kalman filter 
to two dimensions results in excessive computation loads. 
On the other hand, the filters [3-5] are very simple and, 
although it is known that they are not optimal, it was not 
determined yet just how far they are from optimality. There- 
fore, these filters will be compared against the optimum 
non-recursive filter. 

c) The existing recursive image filters exploit only the 
spatial correlation, therefore a three~dimensional recursive 
filter will be developed in order to take advantage of the 
correlation in time. 

d) The possibility of improving the performance of the 
sub-optimum recursive filters by using the non-causal (with 
respect to the direction of recursion) observations closest 
to the estimated pixel (picture element) will be investi- 
gated. The result is a hybrid filter in the sense of using 
some observations recursively and others non-recursively. 

e) An optimum decision rule will be derived for purposes 
of background suppression and target enhancement and subse- 
quent threshold detection. 

f) A conventional Kalman filter will be constructed to 


rack the target centroid from trame-to-frame. 


EZ OVERVIEW 
The image modeling problem is discussed in Chapter II. 
First, a tutorial discussion of the mathematical problem of 


finding the dynamic model of random fields is presented. 


EI 





Second, the results of several experiments with real life 
pictures are presented. Third, a model for a picture and for 
pictures sequenced in time is introduced. 

In Chapter III the sub-optimum two-dimensional recursive 
filters [3-5] are analyzed. The exact computation of the 
error variance is particularly important to those filters, 
because they either do not compute the error variance in 
the calculation of the gains, or use an approximation. An 
exact method is developed to compute the error variance of 
those filters. Using such a method, the filters are compared 
among themselves and against the optimum non-recursive 
interpolator, constrained to the same data set. A recurs- 
ive filter is also introduced that is essentially the 
same as [5], but the computation of gains is accomplished 
without approximation. 

In Chapter IV a three-dimensional recursive filter is 
developed. This filter estimates the pixel gray level of 
pictures sequenced in time by using recursively the obser- 
vations of the estimated frame, as well as those in the pre- 
vious frames. It is an extension of the two-dimensional 
recursive filter [3]. Numerical results are presented to 
evaluate the improvement resulting from the exploitation of 
the correlation in time. 

Chapter V introduces a new class of image filters, 


called hybrid filters. These filters are smoothers that 
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combine optimally the estimate of the recursive filters 
(two or three dimensions) with an arbitrary set of "future" 
observations. Theoretical comparison between the hybrid 
filter and the recursive [3] and non-recursive [6] filters 
is accomplished. Experiments with real life pictures 
are also presented for the purpose of comparison of these 
filters. 

In Chapter VI an optimum decision rule is developed to 
detect targets immersed in cluttered background más 
The target is considered as another texture statistically 
distinct from the background texture. The decision is 
made pixel by pixel and, therefore, it is independent of 
the target shape, but it can also be applied to targets with 
known shapes. The decision is based on the observation of 
the pixel gray level and the background prediction for the 
pixel. This prediction is given by the recursive or hybrid 
filters and may also be given by a non-recursive filter. 
Some special cases are worked out and its performance 
evaluated. The image is modeled as a weighted addition of 
three components: target, background and observation noise. 

In Chapter VII a conventional Kalman filter is constructed 
to track the target centroid (or other points) from frame 
to frame. The target dynamics in the picture is modeled. The 
detector developed in Chapter VI feeds the tracking filter 
With the observations of the centroid spatial coordinates. 


The observation error is analyzed. 


T3 





Chapter VIII presents several results of the simulation 
of a complete target detection and tracking problem, using 
computer generated images. Both recursive and hybrid 
filters are used and compared. 

The final chapter summarizes the results of this 


investigation and presents the conclusions and suggestions 


for further study. 
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II. IMAGE MODELING 


ee INTRODUCTION 

Our objective in this chapter is to obtain a statisti- 
cal model for pictures. Knowing part of a picture, one 
can generally draw certain inferences about the remainder; 
or, knowing a sequence of frames, one can, on the average, 
make a good guess or prediction about the next frame. From 
a Statistical viewpoint, similarity between adjacent pixels 
(picture elements), or , frame-to-frame similarity represent 
a high level of intraframe or interframe correlation. 

Experimental evicence [1]-[2] indicates that a mono- 
chromatic image can be modeled by specifying its value (gray 
level x(m,n) at each spatial coordinate (m,n). An ensemble 
of such images can be modeled by interpreting x(m,n) as a 
random field. 

In this chapter, first, we will address the mathematical 
problem of finding the dynamic model of random fields, given 
the autocorrelation function. Second, we will introduce a 
model for a picture and for pictures sequenced in time. 
Third, we will present experimental results in order to 


validate the proposed model. 


B. DYNAMIC MODEL 
Assume the mean and autocorrelation function of a 
homogeneous (wide sense stationary) random field are given. 


Since we are assuming knowledge of the mean, for convenience, 


ÉS 





we assume the random field has zero mean, therefore auto- 
correlation and autocovariance are the same. 

A highly desirable characteristic of any dynamic model 
is to have it excited by uncorrelated noise. The impulse 
response method always has such characteristics. The diffi- 
culty with this method is that the power spectrum density 
has to be factored, and there is no mathematical theorem 
for factorization in the multi-dimensional case. Since 
there exists considerable evidence, in refs. [1]-[2] and 
in this research, that images are well modeled by auto- 
correlation functions with separable kernels, the difficulty 
of factorization will not arise and, therefore, the impulse 
response method will be used. 

The impulse response method is shown in figure 2.1. 
Assume there exists a system transfer function H(z), such 
that, when driven by white noise W(n), the resultant output 
X(n) 1S a random process with the desired autocorrelation 
function R(k). Using input-output relationships between 
power spectrum densities (PSD): 

S (z) = S,(z)H(z)H(z |) 
O aL 

Since the input is white noise, its PSD is constant, 
S. (z) = K, thus: 


ea 


il 
Se 
S 
II 

pa 


Z2(R(K)) 
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Pigure 2.1 


Filter Response Method 


Therefore, to compute H(z) we have to factor the 
right hand side of the previous equation. This is always 
possible, in the one dimensional case, provided S, (2) is 
a ratio of polynomials. 

In the remainder of this section we present some exam- 
ples to illustrate the method, as well as for later use. 


Example 1 


Consider the autocorrelation function: 


a 
R(k) = ol 2 ok | 


let 


Take the two-sided Z-transform of R(k), 
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Z(R(k)) = R(k)z 
k=-» 
= lose eoa pz+1+ ae er) 
2 2 
Z(R(k)) = — "2 — (2.2) 
(Geez) (loz) 
2 2 
S,(z)H(z) H(z”) A 
E Sa) 


To have a stable model the poles of H(z) must be inside 


the unit circle, therefore we choose: 





l A(z) 

Ez) = E = 
eae l W(z) 

SS er ) 


Thus, the dynamic model is: 
, X (n) = oX(n-1l) + W(n) (223) 


mi tializing equation (2.3) with X(0), and exciting 


2 
), a random process 


with white noise, having variance ar 
X(n) is generated having the desired autocorrelation 


Function. 
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Example 2 


Consider the autocorrelation function: 


2 


g 2 KFL 2 Ik|+1 
R(k) = nn [(1-p, ) ei leer RD ] 
(04 P>) (L+p, 05) DE 1 2 
24) 
where: 
-Q = 
= e 1 r Po = @ 2 pS OZ 


Using results of example 1, the Z-transform of R(k) is: 


Z(R(K)) = A 
(1-9, 2 ) (1-952 ) (1-9, 2) (1-952) 
225) 
where: 
Glide) 
= iba a 2 
Q 140,07 (L Oy la P> ) 


Let's choose: 





Thus, the dynamic model is: 
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X(n) = (p,+p,)X(n-1) — 0,0,X(n-2) + W(n) (2.6) 


It can be easily shown that a particular case of this 


model, when 04 approaches Por 1s: 


R(k) = of of Klis ikl) (2.7) 
where 
2 
3 = 52, o = e! 
1 + 0 


and the dynamic model is: 


ER A A ni (2.8) 


From equation (2.7) we can see that this autocorrela- 
tion function has zero derivative at the origin (k = 0) 

End an inflection point at k = 1/a, provided that p > 0.7 
RR B = a. 

The model given by equation (2.4) is, therefore, a more 
general case and includes, as particular cases, the models 
given by equations (Mia and (2.7). In Elgure 2.2 some 
curves of equation (2.4) are shown for the same "correlation 
time" (the point where the correlation is 37% of the 


maximum). 
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AUTOCORRELATION R(k) 


© 
Oo 


© 
O 


EEN REN 


0.4 





correlation time = 10 


(1) PD, = 0 
(2) P = P3 
2 = 
(8) pP] 20, 
1 
4 
3 
à 
q<_ _A A 4 A 
5 10 15 20 K 
SHIFT 
Eigure 22 


Autocorrelation Functions 
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Example 3 


Consider a two-dimensional autocorrelation function, 
separable in the two dimensions: 
3 | 
Ph (2.9) 


R(17)) er O 


where 


To find the 2-transform we take advantage of separa- 


bility and use equation (2.2) from example 1. 


= 5º (1-0,2) (1-0, 2) 
NS A 
(I-p 2] ) (1-9, 25 ) (1-p_2,) (1-0,25) 
(2.270) 
Let's choose: 
EA = — 7 STO e e 
ROUE ) (1-0, 2, ) 12 72 
Se re + 2 
S; (2,725) = g (l Py Ge Oy, ) 
Thus, the dynamic model is: 
X(m,n) = o, X(m-1,n)+p,X(m,n-1)-p,,p0,X(m-1,n-1)+W(m,n) 
(22.1) 
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Example 4 


Similar to the previous example, consider a separable 


autocorrelation with the kernels of example 2: 


ee > | | 
R(i,JjJ) = o K,(1) K, 0) 622412) 
where: 
1 2 lil+ı 2 jil+l 
K (i) = ut (1l-p Ip = eb ) ] 
y (Piy Poy) CIFP] Pay) 2v ly lv 2v 
RD) = Same equation with i > j 
ven 
The Z-transform is given by: 
HEURES pe azn ey zi 
]=z-» jJ=-» 
= of JC J K Je (kz + 
eae o v l 
| =— C0 J=-» 
= of 2, (K (i)) 25(K, (4) 
pa 2 


Using equation (2.5) it can be easily verified that the 


transfer function is: 
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EZ. 725) 
AZ -l,, mL - | -1 
(1-9) 2) ) (1 09121 pace: P1n29 ye ae 0712) ) 
(MES) 
S. (z,,z.) = 5° OTO 
TA eerie vi 
where: 
Q = 7 PivPav! (y a 
vV i ot Py yP oy T7 2V 
O, = same equation with v > h 
Thus, the dynamic model is: 
X(m,n) = (9,,,+P5,,) X(m-1,n) = Py 99, %(m=2,n) 
+ (944,40 5,) X(m,n-1) = Py pP opt (me n-2) 
+ P40 5,, (0p, +P5,) X(m-2,n-1) 
+ PynPontPintPoy)X(m-l,n-2) 
14722219292 (M-2,n-2) + Wm,n) DA) 


Observe that example 3 is a particular case of this, 


where Poy = = 0. 


Pan 
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To simplify notation, from now on, we will adopt the 


following convention: 


X(m,n) = 


| > 


X + W(m,n) Bel) 


where: 


| 
Il 


row vector of coefficients 


1x 
Il 


column vector of adjacents 


In the case of example 3: 


PE CR pp.) 
X (m~1 ,n) 

x = X(m,n-1) 
X (m-1,n-1) 
E d 


Example 5 


Consider a three-dimensional separable autocorrelation 


function having kernels like examples 1 and 2: 


R(i,j,k) = q? EO RRM a rae) 


where: 
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K, (i,j) = K, (1) K (5) (see example 4) 


Following similar procedures as in example 4, it can 


Be verified that the transfer function is: 





nes ,2,/2-) l 

LE =] = E = a 
(1-9), 2, ~) (1-952, 5) (-p,,2, ) (1=95,25 ~) (1-p,2, 0) 
(2,17) 

E o 

S, (24725723) E g o. (1 Pe ) 

where O, and O, are the same as in example 4. 

Thus, the dynamic model is: 

X(m,n,t) = A X + W(m,n,t) (2:13) 


where: 


E 


| os 


is a row vector whose elements are the coefficients 
of eq. 2.14, followed by these same coefficients 
multiplied by Per and the final element is PE 
is a column vector whose elements are the adjacent 
pixels of eq. 2.14, followed by these same pixels 


located on the previous frame, and the final element 


is X(m,n,t-1). 


In figure 2.3 the adjacent pixels are shown. 
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n-2 nl n n-2 n-l n 


i TIME = t-1 i TIME = t 


Figure 2.3 


Adjacent pixels of example 5 


C. AN AUTO-REGRESSIVE MODEL 

In this section we present a stochastic model for 
images. The basic assumption is that images constitute a 
homogeneous random field with zero-mean (or known mean) 
and known autocorrelation function, separable ın the inde- 
pendent dimensions. These assumptions will be validated 
by experiments in the next section. 

ir Autocorrelation Function 

The autocorrelation function chosen is intended to be 

general enough to include most of the models used in other 
research [2] - [9], as well as to best fit the experimental 
functions measured in this research and ([1]-[2]. The 
hypothesis of separability allows us to examine the auto- 


correlation kernel by kernel. Assume that the process is 


27 





modeled by a kernel of first order, as in Example l. 
The dynamic model is 


X(n) = X (n=l) + W (n) 


ai 
If this model is valid, the modeling error W, (n) must be 
uncorrelated noise. This can be easily verified by computing 
the autocorrelation of W, (n) as follows. We compute W, (n) 


at each point by 


and then compute the autocorrelation of this sequence. 
Computing the autocorrelation of W, (n) is a very good 
test. Assume that W (n) turns out to be correlated. 

In this case we can model W (n) again by a first order 


kernel: 
W (n) = PA + W, (n) 


Assume that, after repeating the same measurements for 
W, (n), we conclude that W, (n) is uncorrelated noise. 
In figure 2.4 it is shown, in the Z-domain, the operations 
that we have performed. 


Thus, the overall transfer function is: 
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Figure 2.4 


Second-order model 


KZ) 
ae = Fa TS 
W, (z) (L-p,277) (1-p,2 


The second-order dynamic model is: 
X(n) = (9, +95) X(n-1) = P1PoX(n-2) T Wo tn) 219) 


We have now a model driven by white noise, which is a 
second-order difference equation. 
The process generated by this model has an autocorrela- 


tion function given by equation (2.4), as we will demonstrate. 
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Z(R(K)) = Q, H(z) H(z 1) 
RE 2) Mo zZ) oz) (1-02) 
1 1 2 2 
AQ, | BO, 
(1-0,27*) (l-p,2) — (l-0,2 >) (1-0,2) 
where: 
p 
ji 
A = EROE a EA 
Pira | 
B = Ei 
(0,797) (1-0,045) 
since 
A = ) = => „IK 
ee Er 1-p 
_ IK | B | k | 
RK = Qi 702 
on PD 
Using R(0) = ge (variance of X(n)) 
Cimo) 
FR 2. Pe id > 
DCE 7149/02) foi So) 
R(k) E ee ee 
GEBET) oo re) Do Ze > 


(2221) 
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Assume, now, that W, (n) is not uncorrelated 
enough. In such a case we proceed to a third-order model. 
Following similar procedures as for the second-order model, 
we have: 


Dynamic Model: 


X(n) = (p,+p,+p,)X(n-1) - (p,p,+0,p3+p,p3)X(n-2) 


+ 0 X(n-3) + W (n) 222) 


ees 


Autoceorrelation function: 











| k | |k | | k | 
oo ee as Dr) 
1 2 3 
where: 
Es 2 z E AZ Re 
A, = 01 (0, 03)! 0793) (1 Pa L P3 ) 
ae ee > = Bee E 
A, = O) (oy 03) (1 0,93) (1 O all P3 ) 
= 2 x aa ee BRENZ 
Variance of W,(n): 
2 
OTe A A ee 
io a 3 
l- £ l- E l- 
21 PD P3 


Sal: 








The general expressions for a k-order model are: 
Dynamic model: 
X(n) = (p,+p,+...pX)X(n-1) =- (01,0540, 034. + +P, _ 1 PR) X(n-2) 


+ (21 PoP3t- +P, _5P)_ 7, ) X(n-3) 


k+1 


+ (-1) 0.7": 


- po, X(n-k) + w (n) (2.25) 


Aueocorrelation function: 











A e ao É 
l 2 K 
Varlance of W (n): 
2 
Q A A 0 (22-27) 
k A A A ° 
l 2 k 
+ +... + 

ee Mo 
a E Pk 


The expressions for the A's are easily inferred 
Brom equation (2.23). 

In example 3 an auto-regressive model is presented 
for a picture with an autocorrelation whose kernels 
in the two spatial dimensions are the first order versions 
(k=1) of equation (2.26). In fact, that model is used in 
most research [2] - [9] with images. In example 4 
an image model is presented whose kernels are the second- 


Order versions (k=2) of equation (2.26). The order of the 
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kernels doesn't have to be the same. Depending on the 
particular picture to be modeled, the order is chosen for 
each direction. 

In example 5 an auto-regressive model is presented 
for pictures sequenced in time. In that case we chose 
second-order kernels in the spatial directions and first- 
order in the time direction. 

At this point, some special features of the proposed 
autocorrelation functions given by equation (2.26) should 
be emphasized: 

a) It includes as a particular case the first 
order model which is used in most research [2] - [9]. 

b) It fits much better many situations where the 
first order model iS a poor approximation. 

c) The resultant dynamic model is auto-regressive 
and driven by white noise. 

d) The cascade feature of whitening the modeling 
error 1S quite simple and adequate for on line parameter 
identification. This feature will be exploited in the next 
item. 

2. Parameter Identification 

Given that the autocorrelation function of an 
image, or a time-frame of images, is assumed to have 
separable kernels as in equatıon (2.26), the next problem 
is to find its parameters. Let's discuss the modeling of 


just a picture, since the extension for the time-frame is 
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straightforward. Having a specific picture, we want to 
identify the parameters of the model. To compute Statistics 
from just one realization of the process we must assume 
the process is ergodic, therefore we consider the picture 
as a typical realization and the averaging in space is 
equivalent to statistical averaging. Also the size of the 
picture has to be sufficiently large. 
Those conditions are quite Severe, and we realize 
that the process has to be restricted to a particular 
class of pictures with common properties. Another way 
of looking at this is to visualize the picture as an array, 
for example, of 100x100, therefore with 10,000 pixels. 
Each pixel can have, say, 100 distinguishable brightness 
values. The number of different pictures is Ore a 
With such enormous numbers of pictures, even if we were 
able to compute the mean and autocorrelation, such moments 
would not be enough information. It would be necessary 
to know statistics of enormously higher orders. The atten- 
tion, therefore, has to be focused on local redundancy, 
encompassing only a few pictures closely related. To have 
a model as accurate as possible, based only in the first 
and second moments, we will need some degree of adaptation. 
Our method is based on the cascade feature of 
equation (2.26). Let's explain the method in detail. 
Assume a picture is given with negligible measurement 


noise. We proceed in the following steps. 
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seep 1 —- Compute the mean and remove it from the picture. 
ERSpP 2 —- Compute the autocorrelation in the horizontal 
direction, with a number of displacements equal to 20% of 
the number of columns. There is no use to go much further 
than this, because the accuracy in the averaging is reduced. 
Repeat the same procedure in the vertical direction. This 
step can be accomplished using the Fast Fourier Transform 
algorithm (FFT). 

step 3 - Having the measured autocorrelation in the hori- 
zontal and vertical directions, assume that equation (2.26) 
is adequate and gives a nice fit. Start with a first-order 
kernel and compute Ply and Pih’ Such that the mean-squared 
error in the fitting iS minimized. The result is a model 


like example 3: 


X(m,n) = Ar Br W, (m,n) 
where: 
le e ne 
x! = [X(m-l,n) X(m,n-1) X(m-1,n-1)] 


Step 1 - At this point we have a model for the picture that 
is first-order in both spatial directions. To test the 
goodness of this model the obvious figure of merit is the 


modeling error W m,n). If the model is adequate, W, (m,n) 


35 








should be uncorrelated noise. To test this, compute the 
autocorrelation of W, (m,n) as in step 2. If it is corre- 
lated in both directions, it means the first-order model is 
not good enough and we have to proceed to a second-order model. 
It could happen that we are satisfied with one direction, 
but not with the other, in this case we have to proceed’ 

to a second-order only in such directions. Assume here 
that W (m,n) is correlated in both directions. 

mtep 5 - Repeat step 3 for W, (m,n) in order to find Day 


and Pop: The model for W (m,n) 1S: 


W (m,n) = A, W, + Wy (m,n) 
where: 
oe oR no 
T 
W = [W, (m-l,n) W., (m,n-1) W, (m-1,n-1)] 


ik i 


The second-order model for the picture is given by 
equation (2.14) in example 4. 
Step 6 - Repeat step 4 for W, (m,n). If it is sufficiently 
uncorrelated, the procedure stops and we have identified 
the order of the model as well as the parameters. Other- 


wise we proceed seeking a higher order model. 
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NOTE: Step 4 is responsible for the decision to proceed 
Se mot and also in which directions, whether horizontal, 
vertical or both. The decision to proceed or not can be 
implemented only based on the variance of W(m,n), by observing 
its reduction after the second-order modeling. If such reduc- 
tionis null or negligible, it means the original model 
is good and it is not necessary to increase the order. The 
decision about the direction to be improved can be imple- 
mented by observing the autocorrelation of W(m,n) at a few 
points, say, the first 5 displacements. 

The computation of the coefficients is quite simple 
because we have reduced the problem to fitting exponentials 
at each step. Such computation is accomplished as follows. 


Assume the measured autocorrelation in one direction is 


even by Y(n), n = 0,1,...K and it is normalized such that 
Y(0) = 1. We want to fit an exponential minimizing the 
mean-squared error, R{n) = e = a Taking the natural 


logarithm of R(n) and Y(n), we want to minimize: 


K 
E” = pino o) = In Y{n)]? 
n=1 
K 
e = ) [-an - In A 
n=1 
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Take the derivative with respect to a and equate 


to zero, the result being: 


K K 
) n In Y(n) -6 ) n in Y(n) 
É n=1 _ n=1 
K K(K+I) (2K41) 
2 
n 
=} 
and 
-q 
p = e`. 


In figure 2.5 a flowchart to implement the proposed 
method of modeling is shown. 

To model pictures sequenced in time the procedure 
is quite Similar. In this case we have a sequence of pic- 
tures, sequenced in time, and apply the method adding the 
extra dimension (time). Example 5 presents a model that is 
first-order in time and second-order in both spatial 
dimensions. 

3. Modeling Of Noisy Images 

In the last item, a method for picture modeling was 
presented in which the picture was assumed to be noise free. 
Such an assumption restricts considerably real life appli- 
cations where there is no knowledge of the original picture. 
This restriction can be removed if the autocorrelation 


Eunction of the measurement noise is known. 
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START 


COMPUTE MEAN AND 
VARIANCE OF THE PICTURE 
REMOVE THE MEAN 


COMPUTE HORIZONTAL AND 
VERTICAL CORRELATIONS 


COMPUTE: 
jy and on 


COMPUTE: w. (m,n) 
and its mean and 
variance. Remove mean. 


COMPUTE VERTICAL 
CORRELATION 





COMPUTE: w; (m,n) 
and its mean and 
variance. Remove mean. 





Figure 2.5 


Parameter Identification Flowchart 
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In what follows, we present the steps to be modi- 
fied in the previous modeling method, in order to take into 
account additive zero-mean white measurement noise, with 
known variance. 


The measurement equation 1s: 
y(m,n) = X(m,n) + v(m,n) 


where v(m,n) 1S zero-mean white noise, uncorrelated with 


X(m,n) and with known variance: 


RECp 2 - Compute the autocorrelation of the noisy picture 
-y(m,n) in both directions, as before, but remove the 


variance R from R.,.,(0,0) to obtain Ry (0,0), Since: 


R,,(0,0) - x (i,j) = (0,0) 
Be (a) 
RR (1,3) # (0,0) 


In figure 2.6 the relation between the autocorrela- 
Elons of X(m,n) and y(m,n) is illustrated in one direction. 
meee à - Since X(m,n) is not known, W (m,n) can't be calcu- 


lated, but its autocorrelation can be computed indirectly. 
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R (0,3 
yy ` J) 


Se) 


j 
Figure 2.6 
Autocorrelation of y(m,n) 
X(m,n) = A] X + W, (m,n) (2,28) 
y(m,n) = X(m,n) + v (m,n) (229) 


Substituting (2.28) into (2.29): 


y(m,n) = A, Y + vim,n) + W (mn) - A 


where: 


K 
h 


[y(m-1,n) y(m,n-1) y(m-1,n-1)|] 


< 
li 


(v(m-l,n) v(m,n-1) v(m-1,n-1) ] 
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mee Ss call: 


vj (m,n) = v(m,n) - A, V (2230) 
u, (m,n) = W (m,n) + v, (m,n) 
y(m,n) = A, Y + u, (m,n) (Za) 


Since y(m,n) and A, are known, u, (m,n) can be 


1 
computed, and, therefore, its autocorrelation in both 
directions. 


Let's examine the autocorrelation of u, (m,n) im 


the horizontal direction. 


Elu, (m,n)u (m,n+j)] = Elv, (m,n)v, (m,n+j)] + E[W, (m,n)W, (m,n+j)] 


l 


E{v, (m,n) v, (m,n+j) J E{v(m,n)v(m,nt+j)] - A,Elv(m,n)v(m,n+}) | 


- A, Elv(m,n+j)v (m,n) ] 


SE 


+ A Elv(men)v (mn+j) JA, 


al 


Ea, nr =, 280,7) 
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0 
E[v(m,n)v(m,n+j)] = |r6(0,j-1) 


0 
S J 
0 
E [v(m,n+J)v(m,n)] = r6 (0,j+1) 
0 
L J 
Elv(m,n)y* (m,n+j)] = | r8(0,3) q r6(0,j-1) 
0 T6(0,]) 0 
r6(0,j+l) O EU, 4 
L J 
Therefore: 
2 2 2 2 l 
= x 
Y A (Ito); +Pih Piy Pih 20,7) 
ode) 
Lh lv rd lh lv 
R (0,3) e OI) T RI W (0740 
SC “i 1 Te 
Br 
2 2 2 2 O. 
a o Sa nio A ne En Jen] = 70 
2 
(0,j) = R (0,1) =~ rp., (1+o O ies 
win, aan Lh Lv 
R4 = Or RE oferty elsewhere 
1 


> 
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In figure 2.7 the relation between the auto- 


correlations of u, (m,n) and W, (m,n) is illustrated. 


R (0,9) 
Sn 
Ea N, 
= | 
-] i 
J 
y + 
Figure 2.7 


Autocorrelation of u (m,n) 


The procedure above is easily extended to calculate 
the autocorrelation of Wo (Mn), as well as higher order 


modeling errors. 


W (m,n) = As W + Wo (m,n) (2252) 
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u; (m,n) = W, (m,n) + V (m,n) 2033) 


Substelietuüte 222) into (2.33): 


u, (m,n) = A, Y + us (m,n) 02.34) 
where: 
u, (m,n) = W, (m,n) + v, (m,n) 12.39) 
vV, (m,n) = vı m,n) - A, Y, 
From equation (2.34) we can calculate u, (m,n), 
since u, (m,n) was computed by equation (2.31), and the 


coefficients A, are given by the fitting of Ry (ey a ce 
To continue the modeling, the autocorrelation 


Of W. (m,n) is derived from R (1,3) uSing equation (2.35). 
Be a 


De EXPERIMENTAL RESULTS 

In this section we present some relevant results obtained 
with real life pictures. Several pictures were analyzed 
and divided in two broad classes: 

Picture A - represents a class of pictures with few 
details (i.e. a low content of high spatial frequency 
Structure). In figure 2.8 a sample of such pictures in 
a three-dimensional plot is shown, where the height (z-axis) 
is the gray level. This kind of plot is very useful to 


emphasize small details. 
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Picture B - represents a class of pictures with many 
details. Figure 2.14 shows a sample of such pictures. 
1. Autocorrelation Function 

The proposed model assumes that pictures have 
Separable autocorrelation functions whose kernels are given 
by equation (2.26). To validate this assumption, the 
autocorrelation of pictures A and B were measured and 
compared with that of the proposed model. Figure 2.9 is 
the measured autocorrelation of picture A, and figure 2.15 
is the same for picture B. Using minimum mean-squared 
error criteria, a first-order model was fitted for both 
pictures, as shown in figures 2.10 and 2.16. Such a model 
didn't work well for picture A, as can be Seen, comparing 
figures 2.9 and 2.10. For picture B it worked quite well, 
as can be Seen comparing figures 2.15 and 2.16. Figure zu] 
is the second-order model for picture A; observe that it is 
a very good fit. Based on these results, we conclude: 

a) The proposed autocorrelation is adequate for 
picture modeling. 

b) The first-order model used in other research 
does not fit well pictures of class A. 

c) Pictures with few details are best fitted by 
a second-order model. 

d) Pictures with many details are adequately 


Fitted by first-order model. 
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2. Parameter. Identification 

The proposed method of parameter identification 
was applied to both pictures. Figures 2.12 and 2.17 are the 
autocorrelation of W (m,n) for pictures A and B, respectively. 
The autocorrelation of the modeling error W, (m,n) is the 
measurement of the "goodness" of the first-order model. As 
predicted before, picture A has a poor first-order model, since 
W, (m,n) is quite correlated, instead of white noise, as 
can be seen in figure 2.12. On the other hand, picture B 
has a nice first-order model, since W, (m,n) is almost white 
noise, aS can be seen in figure 2.17. 

Applying the cascade method of modeling for picture 
A, the second-order modeling error Wo (m,n) was calculated 
and its autocorrelation plotted in figure 2.13. Now the 


modeling error W. (m,n) is quite close to white noise, and 


2 
the second-order model is adequate for picture A, as 
predicted before. 


3. Bandwidth Compression 


Inorder to present another piece of evidence to validate 
the model, it was applied to bandwidth compression. Picture B 


is quantized in 256 gray levels, thus, 8 bits are required 
to transmit the gray level of each pixel. Exploiting 
redundancies in the picture and/or in a time frame of 
pictures (like television), it seems possible to reduce 
the number of bits/pixels to be transmitted, therefore, 


reducing the required bandwidth. 
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Ene model Forlpilcture B as: 


The measured variances of X(m,n) and W, (m,n) are: 


VAR (X) = 49 


VAR (W 


Il 
ER 


1) 


Also, it was seen that W, (m,n) is almost uncorre- 
lated noise. Therefore, instead of transmitting the picture 


X, we can transmit Wy: 


W, (m,n) = 9250) Sek 


Since W is close to white noise, and has much 
smaller variance than X, the quantization levels for W may 
be much smaller than for X. This method of bandwidth 
reduction 1s called DPCM (Differential PulSe-code Modulation). 

Of course, the effectiveness of this method is 
strongly dependent on the "goodness" of the model. 

Figure 2.18 presents the histogram of the modeling 
error W. It can be seen that most "energy" of Wi is 
concentrated around zero. In figure 2.19 the reconstruction 
of picture B is presented, where Wy was transmitted with 


Only 2 bits/pixel (average), instead of the 8 bits/pixel 
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required for X. Therefore, even with a bandwidth reduction 
of 4 times, the reconstructed picture yet carries most of 


the information of the original. 
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Figure 2.8 


Picture A 
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Figure 2.18 


Histogram of First Order Modeling 


Error Of Picture B 
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III. TWO-DIMENSIONAL RECURSIVE FILTERS 


A. INTRODUCTION 

A wide-sense Markov (WSM) sequence, corrupted by addi- 
tive white noise, can be optimally (in the linear least 
Squares sense) estimated by a recursive estimator. The 
fact that image intensity can be modeled by a generaliza- 
tion of a WSM sequence (Chapter II) motivates one to search 
for a two-dimensional Kalman estimator to smooth additive 
white observation noise in images. In references [3], [4] 
and [5] three filters are proposed. These filters use the 
dynamic model in equation (3.2), that is a particular case 
(first order WSM) of the more general model proposed in 
Chapter II. The structure of these filters is also the 
same (see equation 3.4), the difference being the method of 
gain computation. The reason of the existence of three 
Bayesian filters with the same structure, but different 
gains, is that they are all sub-optimum. It has been 
determined [10] that the optimality of the Kalman filter 
is not preserved when generalized to two dimensions, with 
such dynamic models. In reference [9] is proposed 
a vector model, described by a first-order linear n-dimensional 
vector difference equation, that is recursive in one parameter 
and generates the same random field used in [3-5]. With 
such a model the Kalman filter is applied preserving, of 


course, optimality. Unfortunately, the complexity of 
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implementation of such a filter is enormous, due to the 

size of the state vector, that 1s the number of pixels in 

one row, Say, 256 or even 100. Also, its complexity increases 
drastically if we need to extend it to three-dimensional 
random fields, as in the case of time-frame, which will be 
addressed in chapter IV. Similarly, the extension to higher 
order WSM random fields increases its complexity. 

Because of the above considerations, we will examine 
more carefully the sub-optimum filters in [3-5], due to their 
simplicity and adequacy to accommodate an extra dimen- 
sion (time) or a higher order WSM random fields. 

In this chapter we will compare the three filters 
[3-5] among themselves, as well as against the optimum 
non-recursive interpolator, constrained to the same data 
set. We will also introduce a recursive filter that is 
essentially the same as [5], but the computation of gains 
is accomplished without approximation. The figure of 
merit will be the error variance, Since this is the cost 
function used by all three filters. Observe that filters 
[4] and [5] compute error variance in order to calculate 
the gains, but there is an approximation in the recursive 
equation, and this may result in big errors in the gains, 
as well as in the computed error variance. 

In the remainder of this chapter, we will develop an 
algorithm to compute the variance of the estimation error, 
in order to compare the filters having the structure of 


equation (3.4) and dynamic model of equation (3.2). Since 
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filters [4] and [5] have similar algorithms for computation 
of gains, and, because filter [4] is clearly inferior to 


[5], only the former will be used in the comparison. 


B. ALGORITHM FOR COMPUTATION OF ERROR VARIANCE 
Given the dynamic model, the filter structure and the 
gain, we are going to develop an algorithm to compute the 
variance of the estimation error. 
1. Dynamic Model 
The picture is modeled by a zero-mean random field, 


homogeneous, with autocorrelation function: 


== > : 
o uia >, 13 (3.1) 
The dynamic model is: 
X(m,n) = o, 2 (m-1,n)+0,X(m,n-1)-9,,0,X(m-1,n-1)+W(m,n) 
32) 


where the gray level X(m,n) has the autocorrelation of 


equation (3.1) and Wí(m,n) is the random forcing input (or 


modeling error). W(m,n) is zero-mean white noise, uncorre- 
lated with X(p,q), for all pixels (p,q) in region xe 7 
defined in figure 3.1. 

E(W (m,n)) = oW(1-9,*) (1-p,°) = Q (3.3) 
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Figure 3.1 

Definition of the region X. 
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Figure 3.2 


Definition of the region Y 


f 
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2. Filter Structure 


The filters [3,5] use equation: 


> 


em,n) 


where y(m,n) 


where v(m,n) 


[1-K(m,n)] [o X(m-1,n) + ppX(m,n-1) 


= pop, X(m-1,n-1)] t Km ny (m,n) (3.4) 


is the measurement of X(m,n): 


y (m,n) = X(m,n) + v(m,n) BRD) 


is zero-mean additive white measurement noise, 


uncorrelated with X(m,n), and having variance: 


E(v“(m,n)) = Y 


The estimate X(m,n) can be seen as the optimal 


linear combination of the prediction x, mm), and the 


measurement y (m,n): 


mn) 


where: 


Ke > 


(m,n) 


(1-K (m,n) 1X, (m,n) + K(m,n) y(m,n) (2:26) 


o, X(m-=1,n) + o, X(m,n-=1) = 9/92 (m-l,n-1) 


(322) 
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The prediction x (men) carries all the information 
about X(m,n), contained in the "past" measurements. The 
gain K(m,n) is computed such that all the information about 
X(m,n), contained in x, (or the "past" measurements) and 
y(m,n), is effectively used. If that happens, the estima- 


Mion error 
e(m,n) = X(m,n) - X(m,n) 


is uncorrelated with the measurements used to estimate 
em,n). 

It is shown [10] that the estimation error can not 
be orthogonal to all the "past" measurements, but may be 


enly to the region Y defined in Figure 3.2, which 


m-1,n-1' 
excludes row m and column n. 
As a matter of fact, filters [3-5] do not have 


estimation error orthogonal to all measurements in Y 


m-1,n-1 
although incorrectly stated in [3] and [4]. A simple way to 
prove this above is using a counterexample. Let's examine 
the estimation of X(2,3). According to [3] and [4], the 


estimation error of this pixel should be orthogonal to 
Mil,2). Also the error in the estimation of X(2,2) can't 
Memorthogonal to y(1,2), according to [10]. We are going 


to verify that these statements are incompatible. 
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The error in the estimation of x(2,3) 1s 
e(2,3) = X(2,3) - X(2,3) 
Using equations (3.2), (3.4) and the fact that 
W(2,3) and v(2,3) are uncorrelated with y(1,2), it can be 
verified that 
ev.) = [1-K(2,3) [p,e(1,3)y (1,2) 


+ p,e(2,2)y (1,2) = Pyphell,2)y(1,2)] 


Since the filters reduce to the one-dimensional 


Kalman filter for the first row: 


E Spy (1,2) 


ENC Glee ON) 


Thus 


e(2,3)y(1,2) 


Therefore, since in general K(2,3) # 1 and e(2,2) can not 
be orthogonal to y(1,2), we conclude that e(2,3) is not 
orthogonal to y(1,2). 
3. Variance of the Estimation Error 

In the following derivation we will simplify 
notation by dropping the argument (m,n), where this does 
not cause confusion. 

The variance of the estimation error is given by: 

P(m,n) = EI(X-X)?] = E(X?) + E(X?) - 2B (XxX) 
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wets call: 


27 = 2 
B(x) s oe 
E(XX) = b 
P(m,n) = a + g > 2b (328) 


a. Variance of the Estimate 
To compute as observe that X can be written 
as a linear combination of all measurements in region 


Ya a (see figure 3.2). 


f 


X(m,n) = b,,¥(1,1) +... by ¥(1,q) CC Clara) 
b A)... D ; A Coe + 
p1?(P ) Se q) on (p,n) 
BY (m,1) A gy mr) +... Bon Y (Mm, n) 
(3.9) 


Define the column vectors: 


to 
i 
© 
D: 
5 
D: 


|< 
il 


[Y(1,1) y(1,2) ... y(p,1) ... y(m,n)] 
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Eid ci ec ing (3.9): 
x(mn) = Bry (3.10) 
The estimate variance is: 
o? = Elx?) BTE(YYT)B 
or 
9. = B RYB (3.11) 
where: 
Ry = E(X¥") 
Using equations (3.1) and (3.5), it can be 


easily seen that the elements of the autocorrelation matrix 


Ry are given by: 


where i is the sequential number for the pairs 


j for (p,q) as follows 


ee) 


ME) ana 
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(k,2) and (p,q) = (1,1),(1,2),... (m,n) 


Wand 2,22. ma 


p- C Ovar lance- -of X and X 


Using equation (3.10), the covariance b is 


given by: 
_ ee T 
b = E(XX) = B`E(XY) 
er 
a TEE 
b = BR (GFL) 
where 
R = 2 
XY 2 


The elements of the vector Ryy are given by: 
Ryy (1) NN our (3.14) 
where i is the sequential number for the pairs (k,£) with 
(k,2) = (1,1), (1,2),...- (m,n) 


q = Ss mn 


Z1 






meee Calculation of the Coefficients 
Do find the yector of coefficients B, observe in 


equation (3.9) that when: 


a) Sal) 


The result is: 


> 


II 
O” 


Zi mar) a 

It can be easily verified that, using (3.15) 
and (3.4), all coefficients can be computed. By means of 
a Simple computer program, therefore, we can compute the 
exact weight that any measurement has in the estimation of 
K. 

Let's summarize the steps to compute the error 
variance at (m,n): 


r 


(1) To find each coefficient Bo a use equation (3.4) 


^ 


with the values of y given by (3.15). The value of X(m,n) 
will be bo T The initial conditions are: 
f 
X(1,]) = 0 tor 1 = (p-1) or 7 = (q=L) 


(2) Use equations (3.11) and (3.12) to compute EE 
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(3) Use equations (3.13) and (3.14) to compute b. 


(4) Compute error variance by (3.8). 


eee) ANOTHER RECURSIVE FILTER 

In this section we will introduce another sub-optimum 
recursive filter having the same structure as those in [3-5]. 
It is similar to [5], but it computes the gains without 
approximations. 

The objective is to find a way to calculate the gains, 
such that the variance of the estimation error is minimized. 


Let's call the covariance, between X and its prediction 
2 


AN 


Xo’ bo’ and the variance of the prediction E 


Using equation (3.6): 


b = E (XX) = (1-K) E (XX) + KE(Xy) 
But: 
b, = E (XX) 
E(Xy) = E(X%) = q? 
Thus: 
b = (1-K)b + Ko“ le) 


Squaring equation (3.6) and taking the expected 


value: 


13 





E(X*) = (1-K) “E (X ^) Er: 2K (1-K) E (yX) 
Büt: 
A 
A O 
p ( p ) 
Br 
E x = E( XX =p 
(y p’ ( ) D 
Therefore: 
3 = (1-K) “0, 4K" (0 +r) +2K (1-K)b,, (3.17) 


substituting equations (3.16) and (3.17) in (3.8): 


P(m,n) = (1-8) © (0%+0 2D.) + Ker (3218) 


Observe in (3.18) that o“ and >> are independent of 


the gain K(m,n), but they are functions of "past" gains. 


Therefore we can find Rim nm Pas a function OI O a Dy 


p 
Ee and r, such that the error variance is minimized. 
Differentiating (3.18) with respect to K and equating 


to zero: 


of + © é 2b 
K(m,n) = — É er) 


o +o - 2b + Y 
p p 
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The variance of the prediction error is: 


NZ 2 2 
‚ = X-X = T OC 22 
Bean) El ( ) 7] g O b (32.20) 
Ppupstltucing In (3.19): 
E 


Equation (3.20) is intuitively appealing, because it is 
the same as that of the one dimension scalar Kalman filter. 


past tedelng (3.19) into (3:18): 
P = Kr (8222) 


Therefore, with the gain calculated by (3.21), the 
minimum error variance 1S given by (3.22). 

Observe that un is independent of K(m,n), and, 
therefore, can be computed using "past" gains. Its 


calculation is quite similar to that of P(m,n}). 


xX m, = A Y (3:23) 
p` E eS 
where: 
le 
SP Amen 
x RE 2 co violino vim,n-1)]: 
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The variance of Xo is: 


o Zen (3.24) 


where Ry is the autocorrelation matrix of the measurements 


p 
E. 
=p 
The covariance Do mn), Similarly to b(m,n), 1s given 


By: 
b = A Rev (3825) 


where the vector Rey is the correlation between X and E 
De coefficients A are calculated identically to B, 
but using the predictor equation (3.7) instead of (3.4). 
A summary of gain calculation follows: 
(a) Use equations (3.15) and (3.7) to find the 


coefficients Am The initial conditions for (3.7) are: 


Ra =" "O Rep comi Caneca elas: 
(b) Use equation (3.24) to find Ze 
we) Use equation (3.25) to find Do" 
(d) Use equation (3.20) to find Po’ and compute 


m,n) by (3.21). 


D. NON-RECURSIVE FILTER 
The optimum recursive filter must have the same 


performance of the non-recursive Wiener filter, provided 
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the data set is the same. Since two-dimensional recursive 
filters, like those in [3-5], are not optimum, we want to 
find how good they are. 

The minimum-mean-square-error (MMSE) in the estimation 
of X(m,n), using the measurements in region er (see 


figure 3.2), is computed as follows. 


Define the vector 


Ze D) 
y(l,n) h. 
Y = H= j| 5P 
dE AR, As 1 
y (m,n) 
h 
mn 
The estimate of X(m,n) 1s 
X(m,n) = HOY 
and the MSE 
P(m,n) = E([(X-X)*] = E[(X - Hy‘ 
Let: 
E(X°) = oc? 


17 





E 2e T E 
P(mn) = q DER Hie eed (ee 26) 


Differentiate (3.26) and equate to zero: 


U a a) 


Mbstitute (3.27) into (3.26): 


P (M, TL) = ge = HTR 


aan H (3.28) 


AY 
It can be shown that when the observation follows (3.5) 


the MMSE is given by 


E (m,n) = h E SE) 


EE PERFORMANCE OF THE FILTERS 

In this section we present the results of the comparison 
between filters [3], [5] and the one introduced in Section C, 
as well as the optimum interpolator of Section D. 

First of all, we have compared filters [5] and that of 
Section C, since they are essentially the same, but the 
former uses an approximation. This comparison was accom- 
plished by computing the gain and error variance of both 
filters for several situations. The result was that 


they presented practically the same values; therefore 
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we concluded that the approximation made in [5] is very 
good. This conclusion is important, because it allows the 
use of simple recursive equations, developed in [5], for 
gain computation. 

The comparison of filters [3] and [5] against the 
optimum estimator is presented in table 3.1 and figure 3.3. 
In these results the coefficient of correlation is 0.9, 
in both directions, and the variance of X is 1.0. The 
Variance of noise is in the range 0.01-1.0, therefore the 
signal-to-noise ratio goes from 0-20 dB. 

From figure 3.3 and table 3.1, we see that filter [3] 
is better than [5], but both are quite close to the optimum 
filter. In the worst case, the error variance of [3] is 
6.5% greater than the optimum, and [5] is 15%. Also, the 
error variance of [5] is 8% greater than that of filter 
[3]. It was observed, from other results, that [3] and 
[5] approach the optimum for lower coefficients of corre- 
lation and/or low noise. 

Another experiment was performed to change the gain to see 
what MMSE can be achieved with such a filter structure. 
Figure 3.4 and table 3.2 show that filter [3] is quite close 
to the minimum error variance, although it is not exactly 
the minimum as can be verified in figure 3.5 and table 3.3. 
eae conclusion is that filter [3] is better than [5], but 
lt is not the best that can be done. As a matter of fact 
the present method of comparison can always be used to 


compute the best gain, although it is computationally complex. 
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The reason why filter [3] is better than [5] is that, 
although it assumes (incorrectly) that the estimation error 
is uncorrelated with the measurements in region nel 
the gain is chosen such that the error is uncorrelated with 
y(m,n). On the other hand, the error of filter [5] is not 
uncorrelated with y(m,n), although it minimizes the 
up-dated error variance at each point. What happens is that 
the value of such minimum is not only dependent on the 
previous error variances, but also on the covariances of 
the "past" errors, since they are correlated. 

The important conclusion is that the 
M=coursive filters [3], [5] and the one of section C, although 
not optimum, are quite close to optimality. Filter [3] is 
the best and simple enough to be extended to three dimensions, 
in order to exploit recursively the correlation in time. 


This will be done in the next chapter. 
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Noise 
Variance 











Picture Variance 


Picture Correlation 


Table 3.1 


Optimum 


-0083729 
-0262313 
.0477314 
.0706634 
.0940204 
.1172805 
. 1401542 
„10624769 
.1841581 


„2057533 


-0 





= 0.9 
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Filter [3] 


.0084204 
202690593 
.0497128 
.0743120 
.0994702 
.1245305 
- 1491169 
. 1730206 
‚1961201 


.2184254 


Filter [5] 


.0083936 
.0266656 
.0494879 
.0747814 
.1014409 
‘1287582 
1562232 
.1834437 
+2101098 


- 2359806 








ERROR VARIANCE 


es 


PICTURE VARIANCE = 1.0 
PICTURE CORRELATION = 0.9 
(0) - optimum filter 

(*) = filter [3] 

(+) - filter [5] 





NOISE 


(Standard deviation) 


Figure 3.3 


Comparison of the Filters 
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Gain 


.07982749 
1095379301 
11175853 
.12772405 
. 14368957 
213965503 
217/962050 
719138608 
220755148 
„223531706 


„23948258 


Picture Variance 


Noise Varlance 


Picture Correlation 


.0 


83 


Table 3.2 


Error Variance 


0. 
0. 


0. 


Ors 
0% 
0. 
0. 


O. 


.9 


26742077 
24629241 


293272502 


„22441578 
219953461 


21542539 


21918494 
22175616 
22571 TIT 
23097461 


23713821 





Potter [3] 


Falter [5] 











ERROR VARIANCE 


PICTURE VARIANCE I; 
NOISE VARIANCE Je 
PICTURE CORRELATION = 


26 
24 


E aan TT 
0.200 


0 
e 


ponlo 
o 07s 190 125 150 8 175 


Gain 


Figure 3.4 


Error Variance vs Steady-State Gain 
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Gain 


023339393 
0. 30999994 
0319939993 
0732399992 
07329993991 
023493993990 
0.35999990 
036999995 
0.37999994 
02383939933 


099939939932 


Table 3.3 


Error Variance 


0.09970987 Filter [3] 
0.09923345 
0.09892100 
0. 09875959 
0.09873748 Best Gain 
0.09884429 
0.09907085 
0.09940886 
0.09985095 
0.10039049 


0.10102153 Filter [5] 


Picture Variance = 1.00 


Noise Variance 


= 0.25 


Picture Correlation = 0.90 
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Figure 3.5 


Error Variance vs Steady-State Gain 
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IV. THREE-DIMENSIONAL RECURSIVE FILTER 


A. INTRODUCTION 

The two-dimensional recursive filters, presented in 
Chapter III, use only the spatial correlation between 
pixels in order to estimate the gray levels of noisy pic- 
tures. In this chapter we introduce a recursive filter 
that takes advantage of both correlations in space and 
in time, when one has a group of pictures sequenced 
in time. Experimental evidence, presented in [1], shows 
that television pictures are correlated in time. This is 
also an intuitive result, since knowing one frame we often 
can guess the next. 

The filter developed here is an extension of the two- 
dimensional recursive filter [3], since it is the best, 
according to the results of Chapter III. In any case, its 
performance will be analyzed, in order to evalute the 


improvement resultant in using the time correlation. 


er FILTER DESIGN 
The filter will be developed under the following 
Bonditions: 
(a) The time-frame is modeled as an homogeneous random 
field with zero-mean (or known mean) and autocorrelation 
function R(i,j,k), where i,j and k are the distances between 


pixels in the vertical, horizontal and time coordinates, 
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respectively. The autocorrelation is like the one intro- 
duced in Chapter II, with kernels given by equation 
7-26). 

(b) The dynamic model is given by the partial difference 


equation: 


2 (m,n,t) = AX + Wim,n,t) (4.1) 


Mete A is a row vector of coefficients and X is a column 
vector of adjacents. The modeling error W(m,n,t) is white 
noise and uncorrelated with the X's in region X ant (see 
figure 4.1). 

In order to make clear the present derivation, let's 
assume a specific form for the autocorrelation, say, first 
order in all dimensions. The development remains valid for 


the k-order kernels given by equation (2.26). 


The autocorrelation is: 


R(i,j,k) = s? hd Q 3 o, E .. Logus = Oi 
dc 

The coefficients are: 
Bees le, Ph Pe ~Oh TPP, Py Pe PyPyPr} (4.3) 


The adjacents (see figure 4.2) are: 
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Figure 4.1 


Definition of the region X the region 


m,n,t' 
includes the past frames and the pixels shown 


m-L 





Figure 4.2 


Configuration of the adjacents 
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x" = [X(m-1,n,t) Mim, n= 1, €) X(m,n,t-1) “nel ,na-1l,t) 
(4.4) 


The variance of W is 


2 
= ) (1-47) (1-p,.°) 


E [W“(m,n,t)] = Q = o? (1-0, É £ 
(4.5) 


(c) The pictures are contaminated with zero-mean white 


noise, uncorrelated with the X's. 


(4.6) 


.Y(m,n,t) = X(m,n,t) + v(m,n,t) 


The noise variance is: 


Deren) = r 


l. Filter Structure 


Since we want to generate a random field X as close 


as possible to X, let's choose a structure similar to the 


dynamic model (4.1): 


(4.7) 


X(m,n,t) = B x % AE Am, Me) 


In what follows we will drop the argument (m,n,t), 


where it could not cause confusion. 
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Observe that X can be written as a linear combina- 


tion of all measurements within region Ya defined in 


ne 
Figure 4.3. The optimum recursive filter must have the 
gains such that the weights of the measurements are exactly 
the same as for the non-recursive estimator, constrained 
to the same data set Yint’ The criteria of optimality is 
minimization of the Mean Square Error (MSE). It will be 
seen that it is impossible to have an optimum recursive 
filter, constrained to the dynamic model (4.1). 
2. Orthogonality Principle 

A necessary condition of optimality is that the 

estimation error be uncorrelated with the data set. Let's 


apply this condition to find the coefficients B and K. 


The estimation error is: 

e(m,n,t) = X-X (478) 
The orthogonality condition is: 

Ele(m,n,t)y(p,q,r)] = 0 (4.9) 


Bor all (p,q,r) in Y (see figure 4.3) 


NE 
First, let's apply (4.9) to the measurement y(m,n,t): 


E(ey) = E[(X-X)y] = Ely (X-BX-Ky) ] 


E [Xy-Ky“-ByX] 
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time 


Figure 4.3 


Definition of the region ee 
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Using equations (4.2) and (4.6): 


E(ey) = (1-K)0? - Kr - BE(XX) = 0 


K(o2+r) + BE(XX) = 0° (4.10) 


Let's apply (4.9) to the rest of the measurements. 


Amos ertuting (4.7) into (4.8): 


e =X-X=X-Ky-BX 
Z oo S B XI Be 
= 1 
e = (1-K) [X - Tx B X] = KVOT BTE Ca 


substituting (4.11) into (14.9): 


BE Iy(par)el + (1-K)E[(X- B X)y(p,q,r)] = 0 
4,12) 
er esq rr * vim,n,6), andoin ae 
Ene problem now is to mchoose te coefficients B 
and K, such that equations (4.10) and (4.12) are satisfied. 


At this point we follow [3] by choosing 
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With such a choice, the second term of (4.12) 


becomes: 


AENA O, E) Y (Pr) ] 


Where we have used (4.1), and this term vanishes 
because W(m,n,t) is uncorrelated with all measurements in 


region X nt” Substituting (4.13) into (4.10): 


5° + b 

K(m,n,t) = = (4.14) 
o" + b +r 
P 
where we have defined: 
Zement) = PA X 
D a a 
(24383175) 

Bm se = E(XX 
p | n,t) ( p? 


Equation (4.15) seems to be an adequate choice for 
the first-step predictor, also (4.14) is an intuitive result 
for the gain, since it is unity for zero noise (r=0) and 
decreases when noise increases. 

Unfortunately, the first term of (4.12) remains and 
we could not force it to vanish. Choosing K by (4.14) we 
have forced the estimation error to be uncorrelated with 
IO, E), but it is not correct to conclude, by induction, 
Siam cic estimation errors of the adjacents (Vector e) are 


uncorrelated with the other measurements. 


94 





Take, for example, the error: 


N 


e(m,n,t-1) = X(m,n,t-1) - X(m,n,t-1) 


The estimate didn't use any measurements of frame 
t, therefore the error can not be uncorrelated with them. 

As a matter of fact, the errors might only be 
uncorrelated with the measurements within region ee 
In reference [3], for the two-dimensional case, it is 
stated that the error is uncorrelated with all measurements 
in a similar region. Unfortunately, that is not correct. 
The estimation error, in the two or three dimensional cases, 
is uncorrelated only with y(m,n,t), because K was computed 


under such conditions, and with the measurement y(1,1,1). 


This can be seen from (4.12) and (4.13): 


Biley(p,g9,2)] = (i-K)EIA e y(p,q,r) } (4.16) 


tor (4,2) = (a,a, 2) 


Initializing the filter with the mean value of X 
(=0), it can be verified, using (4.16), that e(m,n,t) is 
uncorrelated with y(1,1,1). 

Although equation (4.14) is not the optimum choice 
for the gain, we have shown in Chapter III that it was the 
best in the two-dimensional case, therefore we hope it is 
also good in three dimensions. This will be verified at 


the end of this chapter. 
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Using (4.7) and (4.13), the filter equation becomes: 
alee) eee PAX ok KY (4.17) 
Alternatively, using (4.15): 


(mane). = (1-8) X, + Ky (4.18) 


x > 


eB. Gain Computation 


Let's develop a recursive equation to compute Bo 


and, therefore, the gain. Define a function: 


bar (mnt) = E(X(m,n,t)X(p,q,r) ] 
Using (4.17): 
P oar Esto, gocA E + Kip,d, =) Rim-—p,n-c,e-4r) (4.19) 
where: 
FT = [F F F F F 
= le, E o ely teat e plg, 


Fp,q-l,r-1 *p-1,q-1,r-1! Gao 


^ 


Initializing with X equal to the mean of X, outside 


the picture, the initial conditions for F are: 
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rence ce = 0 For pr or ES O 


Now, using (4.15) and (4.19): 


b, (m,n,t) 


I 
| > 
Ir 


ae) 


where Fo PomEne value OL E for (p,q,r) = (m,n,t). 
Summary of gain calculation: 


(a) Use (4.19) to compute each component of the vector 


(b) Use (4.21) to compute Do tmn t). 
(c) Compute K(m,n,t) by (4.14) 

Observe that the calculation of each gain K(m,n,t) 
requires scanning the time-frame from (1,1,1) through the 
pixel (m,n,t). Fortunately, as will be seen, after a few 
frames the gain reaches steady state, therefore the compu- 
tation can be reduced to the top left corner of a few 


frames. 


C. PERFORMANCE OF THE FILTER 

The method used in Chapter III, for the two dimensional 
case, can be easily extended to three dimensions. The 
Variance of the estimation error is: 


P(m,n,t) EL(x - x)?; 


= E(x?) + E(X2) - 2E(XX) 
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Or.: 


Plm, nae) = q +79 - 2b 4222) 


where: 


B(X*) = q? 
22 = 2 
E(X) = I. 


ti 
< 
= 
I 
D 


The covariance b can be calculated using Do and K: 


Sm, n,t) = E (XX) = (1-K) E (XX) + KE (xy) 


bimm,t) = (1-K)b, + Ko? (4.23) 


To compute the variance of the estimator, write X as 
a linear combination of the measurements in region Y 


nt’ 


de) ce (4.24) 


where -€ is a column vector of coefficients, function of 
fee gains, and Y is a column vector containing all measure- 


ments in Y 
mnt 


The variance of X is: 


IS 





g = i RC (425) 


where Ry Sine correlation matrix of the vector Es 


_ T 
Ry = E(YY”) 
The key, in this computation, is the method of finding 


the coefficients c. This can be easily accomplished 


using the filter equation (4.17) with: 


ee dor oe ee GRR Cio dE) 
yli,j,k) = 
O ta sn (ao See Ses 
Een, = 0 for k = (r-l) or 
1 = (p-1) or 
Io 


Initializing (4.17) as above, and starting at (p,g,r), 
the resultant value at (m,n,t) will be the coefficient of 
y(p,q,r) in the estimate X(m,n,t). In this way all coeffi- 


Clents can be found and, therefore, the error variance. 


We are also interested in the first-step predictor 


Given by (4.15). Its error variance is: 
Bim mn, EMP E(X x ae (4.26) 
p p 
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It can be related with P(m,n,t) uSing (4.18) and (4.26): 


X-X = X- (ITK) X5 - Ky 


LO) CRE x, o 


Squaring and taking expectations: 


P = (1-K) P +kK er 
p 
Po(myn,t) = À (p- xr) (4.27) 
P (1-K) 


RR RESULTS 

In figures (4.4) through (4.9) some results of gain 
calculations are shown. For the first frame the gain is 
the same as in [3] for two dimensions. From these results 
we can observe that the gain reaches a steady-state value 
very fast, at about frame number 4. Also, for the same 
frame, the steady-state is reached at about pixel (4,4). 
That 1S an interesting result, since it simplifies substan- 
tially gain calculations. Also we can think of using a 
constant gain for the whole time-frame, or a few gains 
for the first 4 frames. 

In tables 4.1 through 4.3 a comparison is shown 
between the three-dimensional filter, estimate and prediction, 
and the two-dimensional recursive filter. It can be seen 


in table 4.1 that exploitation of time correlation can be 
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quite advantageous in cases where time correlation is 

very high, compared with spatial correlations. In this 
example, the estimation error variance is reduced by about 
33%, due to the use of time correlation. In table 4.2 

a comparison with equal amounts of correlation in time and 
space is shown. In this example the improvement was about 
10%. It can be seen in table 4.3 that exploitation of time 
correlation is not so advantageous when the images are 

more correlated in space than in time. In this example 


the improvement was only around 6%. 


BOT 
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Figure 4.4 


Gains for the First Frame 
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Gains for the Time-Frame 
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Figure 4.6 


Gains for the First Frame 
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Figure 4.7 


Gains for the Time-Frame 
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Figure 4.8 


Gains for the First Frame 
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Figure 4.9 


Gains for the Time-Frame 
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Table 4.1 


Three-Dimensions Two-Dimensions 


Noise 






mee 2493 
.212 .450 
; : ; 132 397 
Picture Variance = 1.00 
Spatial Correlation = 0.70 
Time Correlation = 0.95 


Table 4.2 


Three-Dimensions Two-Dimenslons 
Estimation | Prediction Estimation Prediction 








Noise 
Variance 







Picture Variance = 1.0 









Spatial Correlation = 0.8 






Time Correlation = 0.8 
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Table 4.3 


280 







Noise 


Variance 





.301 . 427 
250 . 268 - 400 
229 367 
180 .324 | 
«114 203 


—— s 









Picture Variance = 1.0 







Spatial Correlation = 0.8 






Time Correlation = 0.7 
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V BRIO ETLTERS 


A. INTRODUCTION 

In this chapter a new class of image filters is intro- 
duced, called hybrid filters. The recursive filters, 
presented in Chapters III and IV, utilize only part of the 
data set, arbitrarily defined as "past" measurements. Since 
the whole picture is often available, an optimal smoother 
might be sought in order to utilize all the information 
available. For the two-dimensional recursive filters of 
Chapter III, the smoother would be the optimum combination 
of the four estimates obtained by scanning the picture 
starting at each corner. The first difficulty is that 
those filters are not optimum resulting in a sub-optimum 
smoother. Second, the smoother would require scanning the 
picture four times, therefore increasing substantially its 
complexity of implementation. 

The hybrid filter introduced here is a smoother that 
combines optimally the estimate of the recursive filter 
(two or three dimensions) with an arbitrary set of "future" 
measurements. It will be applied in picture enhancement 
and compared against the recursive filter as well as the 
non-recursive filter presented in [6]. Some examples of 
hybrid filters designed to predict the pixel gray level will 
also be presented. These filters will be applied in 


Chapter VI for purposes of target detection. 


IE 





The hybrid filter is particularly useful in some 
applications where the pictures are not from a fixed 
scenario. For example, the camera has a movement in order 
to follow some target. In such cases, the background may 
change so much that past frames don't carry enough informa- 
tion about the present frame. In this case, the three- 
dimensional recursive filter can't be applied, but a hybrid 
filter can be designed using recursively the measurements 
in the present frame, and, non-recursively, the "crude" 


measurements in the previous frame. 


B. TWO-DIMENSIONAL HYBRID FILTERS 

In this section we present two examples of hybrid 
filters in two-dimensions. One is a smoother designed to 
enhance pictures contaminated with additive white noise. 
The other is a filter designed to predict the pixel gray 
level, based on noisy measurements of previously scanned 
pixels. Here we assume the picture is scanned row by row 
Erom top to bottom. 

l. Two-Dimensional Hybrid Smoother 

The recursive filters presented in Chapter III 

utilize the measurements in region Yan (see figure 3.2) 
to have an estimate a To improve this estimate we 
are going to use also a selected set of measurements in 
the neighborhood of (m,n), say Y, and combine optimally 
with un), in order to have a better estimate x, (m,n). 
Let's choose for Y the set of five neighbors shown in 


megure 5.1. 
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n+1 
In n ae 
m-l e ‚08 x 
/ 
m / O , Y¥ 
TA i Fa 
m+1 /o o 0: 
x er A 


Figure 5.1 


Two-dimensional hybrid smoother 


Krqure 5.2 


Non-recursive filter 


12 





Thus 


X, (m,n) = hg tre | 
where 
Xm) Ay 
wos a) A, 
y (m+1,n) h, 
a h = 
v(m ley ral } hy 
y (m,n+1) he 
pme] Ae 


Let's calculate the weights h in order to minimize the 


mean-square-error in the estimate. The error is 
e(myn) = X-X, = X- pl 
r 1 ng 


The error varlance is 


I 


P(m,n) = El(X - h g)“] (5.2) 


PSS rentiating (5-2) with respect to h and equating 


to zero 


2 RTE (5.3) 


10, 3 





where 


Elgg.) 


R 
ga 


Ryg E (Xg) 


The measurement equation is 
y(m,n) = X(m,n) + v(m,n) (5:47 


where v 1S zero-mean white noise, uncorrelated with X, 


having variance 
2 _ 
Elv” (m,n)] E 


The gray level X is a zero-mean (or known mean) 


homogeneous random field with autocorrelation 
R(i,j) = ElX(k,2)X(k+1,2+3)] dB 


The elements of the autocorrelation matrix E that 
don't include an are easily found using equations 
(5.4) and (5.5). To find the elements which include x, 
observe that in Chapter III a method was developed to cal- 
culate error variance of recursive filters [3,5], where 


2 


E(X) is also computed as well as E(XX). It can be 


verified that the same method can be used for higher order 
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recursive filters, though in this derivation we restrict 
ourselves to recursive filter [3]. 


Therefore, using (3.11) and (3.13): 


E (X2) 


$ 
a 


(5.6) 


E (XX) 


ii 
O” 


It can be verified that the covariances between 
X and the measurements Y can be expressed as a function 


of b given by 


Ely (m+1,n=1) ] = Ph? 

PAC E ob 

Era] = P PhP (5.7) 
E(Xy(m,n+l)] = pb 

E dO] = PP nº 


Therefore we have 


IES 





E 
o 


Fe PPRP PP PP RP PhP o 
o tr ze 252 2,2 2 2,2 
Ph Ph PPh Py Ph 
J rr a O 2 Be 
Ph Ph°v Py Ph 
Rg É 
2 2 2_2 
o”+r PI Py O 
2 2 
Oo Fr PS 
o tr 
(5,8) 
b 
52 
PPh 
5° 
Py 
R = 
XG 2 
Py PH? 
57 
Ph 
Ee 
PPh 
where we have made use of 
ae 2 al 
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a. Performance and Comparison 
The error variance of the hybrid smoother is 


calculated substituting (5.3) into (5.2): 


P (m,n) 


N 
a 
l 
fo) 
bd 
u 
O 


At this point it is interesting to compare the 
hybrid smoother and the recursive filter [3], since we 
know how to calculate the error variance for both. It is 
also interesting to include in such comparison the non- 
recursive filter presented in [6]. This filter uses a 
small window of 3x3, shown in figure 5.2, to estimate the 
pixel gray level at the center, minimizing the mean-square- 
error. The hybrid smoother must have superior performance, 
since it utilizes the same measurements of [3] and [6] 
together. However, we can't say which is better whether 
[3] or [6], since they use distinct data sets. Although 
filter [3] uses a larger data set than [6], the former 
uses all the closest adjacents to the estimated pixel 
m,n). 

In table 5.1 some numerical results are shown 
to help the comparison. As expected, the hybrid filter 
presented the best performance with error variance about 
25% smaller (for p = 0.94) than filter [6]. The non-recursive 
filter [6] presented better performance than the recursive 
filter [3]. The error variance of the optimum interpolator 


constrained to the observations in Ya = is also included in 
f 


ET 





Table 5.1 


ce varie sid 
Hybrid Filter [6] Filter [3] 












Optimum 
Interpolator 


Correlation 


0.400 009599 
02357 0359 
0.300 0.294 
0.218 0205 
0197 omgi 
OES 0=1.55 
0.146 0.124 








Picture Variance = 1.0 










Noise Variance = 1.0 





this table. It is interesting to see that the use of just 

5 non-causal observations by the hybrid filter were enough 

to fully compensate the suboptimality of the recursive 

filter [3]. Observe that the 9-point non-recursive filter 
[6] also presented better performance than the optimum 
interpolator, even for spatial correlation as high as 0.94. 
This indicates that most of the information about X(m,n) 
resides in its nearest neighbors. This motivates the deriva- 
tion of simple filters like the hybrid filter introduced 
here. It is much simpler than the optimum recursive filter 


[9] with insignificant loss in performance. It also has 
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the same order of complexity as the non-recursive filter 
[6], but superior performance. Complexity here means the 
number of multiplications, which is 9 for the non-recursive 
Mercer and 10 for the hybrid filter. 

The hybrid filter may be seen as a way of 
artificially increasing the window of filter [6]. 

b. Experimental Results 

In this section we present the results of some 
experiments where the filters were used in picture enhance- 
ment. In this experiment, the picture of figure 5.3 was 
contaminated with white noise. Assuming knowledge of the 
noise variance, the picture was modeled by the method 
presented in Chapter II. With this model and noise vari- 
ance, we designed the three filters: hybrid, recursive 
and non-recursive. The theoretical performance was computed 
for each filter, as in the last item. A signal-to-noise 
ratio (S/N) was defined as the ratio between the variance 
of the original picture (figure 5.3) and the variance of 
the noise. A processing gain was defined as the ratio of 


(S/N) 's at the filters output and input: 


(S/N) | 
F = 
(S/N) ; 
Table 5.2 summarizes the results. It can be 


seen that the computed performance is pretty close to the 
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results of the experiments. The hybrid is clearly superior 
to the others and the non-recursive filter is better than 
the recurSive, except for the very low S/N of -3 dB, when 
the former was a little better. However the subjective 
evaluation of the pictures (figures 5.4 through 5.15) is 
more difficult. In our evaluation the hybrid filter had 
superior performance, therefore validating the theoretical 
results, as well as the measurements of table 5.2. However, 
the results of the non-recursive filter, in our subjective 
evaluation, were inferior to the peewee filter in all 
three experiments, contradicting the measurements of table 
5.2. Probably the mean-square-error is not a good figure 
of merit for visual perception, however there is no doubt 
that considerable improvement was achieved by processing 
the noisy pictures using such criteria. 
2. Two-dimensional Hybrid Predictor 

Assuming that the picture is scanned row by row, 
we may define "past" measurements as those pixels already 
Scanned. Then the gray level of pixel (m,n) may be pre- 
dicted using "past" measurements. In figure 5.16 the 
configuration of a hybrid filter designed to predict X(m,n) 
is shown. This filter combines optimally the prediction 
x (m,n) of the recursive filter and the measurements Y, 


P 
shown in figure 5.16. The result is an improved prediction 


A 


X(m,n) that utilizes all the information contained in the 


measurements Zan and Y. 
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Figure 5.5 


Image of figure 5.4 filtered by the hybrid filter 
(SNR = 11.5 dB) 
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Figure 5.6 


Mage of fiqure 5.4 filtered by the recursive filter 
(SNR = 9.3 dB) 
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Pigure 5.7 


Image of figure 5.4 filtered by the 9-point non-recursive filter 
(SNR=10.4 dB) 
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Image of figure 5.8 filtered by the hybrid filter 


(SNR = 9.8 dB) 
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Image of figure 5.8 filtered by the 9-point non-recursive filter 
(SNR = 8.2 dB) 
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Image of figure 5.12 filtered by the hybrid filter 
(SNR = 8.1 dB) 
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Figure 5.14 


image Of figure 5.12 filtered by the recursive filter 
(SNR = 6.6 dB) 
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Three-Dimensional Hybrid Smoother 
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The design of this filter is similar to the hybrid 
smoother of the previous item, therefore, we present below 
the results. 


The autocorrelation function used is 


ee | | 
R(1,3) = © P3 Ph 52.10) 
The predictor is 
X(m,n) = hig 
where 
f LR 
Xx, (m,n) Ay 
yım= 1, Tl) no 
g = y (m-1,n+2) a na 
vm do mi) hy | 
= | 
mL nr) he | 


The weights h are given by 


= 
h = R R 
== ag Xg 


Using the results of Chapter III we can compute 


AT the elements of R and R. . Let's call: 
gg Xg 
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>; A 
E LX A = oO 
[ p (m,n)] p 
t-il) 
XX = b 
E[ pl p 
It can be verified that the cross-correlation 
between XD and Y can be expressed as a function of Dor 
as follows: 
Elx, y (m-1,n+1)] = PPP py 
EIX v(m-1,n+2)] = pp T 
p ae ie 
(5.72) 
_ 2 
ElX y (m-2,n+1)] > Pe OP, 
Pe aa nna] ns 
p í wen BD 
Using equations (5.10), (5.11) and (5.12), we can 


Calculate the autocorrelation matrix Bea and the correlation 


vector R 
x 


Ze PyPhlp 00 Po o, Pub, Dy Ph Do 
gee 0,0? o 0 o. py 0 
icc = Str oo O 
of+tr on 
of +r (5.13) 


m 


PO 





Xg v°h 


C. THREE-DIMENSIONAL HYBRID FILTERS 
In this section we present two hybrid filters 
applied to pictures sequenced in time. The first filter 
is a smoother to be used in image enhancement and the 
other is a predictor. The design of these filters is 
quite similar to those of the previous section. 
l. Three-dimensional Hybrid Smoother 

The filter structure is shown in figure 5.17. 
The recursive filter is that developed in Chapter IV. The 
estimate ne) of the recursive filter is optimally 
memotned with the set of measurements Y, in order to have 
a better estimate x, (m,n,t). The design is quite the same 
as in the two-dimensional case, thus we present only the 
results. 

The autocorrelation used for the time-frame is 


R(i,j,k) = of 0 (5.14) 
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The estimate is 


A(m,n,;t) = h'g 
where: 
Ve Misty te ye) h, 
ini ne) A, 
y (m+1,n+1,t) hy | 
| 
ym ARE) he 
a y (m-1,n+1,t) h = he 
Van inde 10) h- 
PAIN El) he 
y(m+1,n+1,€-1) | Pi) 
x | | 
PEL, EC 1) | Aig 
| h 


yml ntl tD | 


= 
| 
| 


melweights h are given by 


aa 
R R 
g9 XI 


Using the same notation as in Chapter IV, we have 


E[X*(m,n,t)] = o 2 
e 
es 


E[XX] = b 
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ZN 


As before, the correlation between X and the measure- 
ments Y can be expressed as a function of b as follows 


Mi (m,n,t)y (m+1,n+],t+k)] = op b (216) 


where 


ir ris) = A UA ON O RO), 
ld 00) (SA) OL) / (a) , 


A CL) 


Using equations (5.14), (5.15) and (5.16), we can 
calculate the autocorrelation matrix = and the correlation 
vector a then the weights h. 

2. Three-dimensional Hybrid Predictor 

This filter combines optimally the first-step 
prediction x of the two-dimensional recursive filter [3] 
and the set of measurements Y shown in figure 5.18. Since 
its design is quite similar to the two-dimensional case, 
we present only the results. 


The autocorrelation function is given in (5.14) 


and the hybrid prediction 1s 


Dm, E) 


[ 
[D 
ka 


where: 
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Figure 5.18 


Three-dimensional hybrid predictor 
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The weights h are given by 


Similarly to the two-dimensional case: 


a 


E [X 
p P 


ll 
Q 


Vol) 


ll 
O” 


E [XXS] 


AN 


And the correlation between Xo ande‘ istgiven Dy 
RR) ms ny, 21) = p ES 13] Do (ee) 


Vv Ph Pe 


where 
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(1,9) = ee eon O, =1),(07,0),10,L) +, 


A , (1,0) , (1,1) 


The autocorrelation matrix Rag and the correlation 
vector R, can be calculated using equations (5.14), (5.17) 
g 


and (5.18), then the weights h are computed. 


D: COMPARISON OF THE PREDICTORS 

The predictors developed here will be used in the target 
detection problem addressed in Chapter VI. In this section 
we present some numerical results in order to evaluate its 
performance. 

In Tables 5.3 through 5.5 the error variance of the 
recursive and hybrid predictors for several situations is 
Shown. From these results we can see that the exploitation 
of the time correlation can be quite advantageous, mainly 
mor the case of high correlation in time (table 5.3), where 
the error variance of the three-dimensional filters is around 
30% of that given by the two-dimensional filters. We also can 
observe in the three tables that the 3D-hybrid filter is 
better than the 3D-recursive filter for observation noise 
Variances above 0.2, although the former exploit only the 
time correlation with the previous frame. This result rein- 
forces the basic idea of the hybrid filters, that is the use 
of all the neighbors of the pixel (m,n), including those 


that are not causal. 
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Table 5.3 


Error Variance 


three-dimensional two- | two-dimensional | 
hybrid recursive o. 


2902 












Noise 












Variance 





Table 5.4 


Error Variance 


Noise three-dimensional two-dimensional 


Jol 303 427 364 


.400 .340 
00 EZ 
24 2710 
203 227 


Picture variance = 1.0 
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Table 5.5 





Error Variance 













Noise three-dimensional two-dimensional 
-3/6 53308 -427 - 364 
.400 .340 






Picture variance = 1.0 
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VI TARGET DETECTION 


A. INTRODUCTION 

In Chapter II we have developed models for images and 
for a sequence of images in time. Using such models we 
have constructed two and three dimensional recursive and 
hybrid filters to smooth out noise in images. In this 
chapter we are interested in the detection of targets using 
low contrast images. This problem is usually referred to 
as background suppression, where background means a kind 
of texture that predominates in the picture. In this con- 
text, the word target means the information that we are 
interested in extracting from the pictures. It may be a 
ship, an airplane or some kind of texture statistically 
different from that of the background. 

Both textures, target and background, are modeled as 
in Chapter II, but it is assumed that less "a priori” 
knowledge is known about the target. The basic idea in 
Beeection is to use the prediction of the image filters 
previously developed. Knowing this prediction and the actual 
observation of the pixel we develop here, a decision rule 


to decide to which texture the pixel belongs. 
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B. TARGET AND BACKGROUND REGIONS 

Figure 6.1 shows two distinct regions in the picture: 
background and target regions. In the visible spectrum 
what happens is a replacement process, we have background 
or target. On the other hand, in the infrared spectrum 
what seems to happen is some combination of the radiations 
from the target (e.g., a plane) and the background(e.g., 
clouds). To include these situations, let's assume that the 


gray level in the target region is given by 


fam) = ax(m,n) + T(m,n) + v(m,n) (6al) 


This model assumes alinear combination of target and 
background, normalized with respect to the target, and also 
consider an additive white observation noise v(m,n), 
uncorrelated with both target and background textures. The 
background texture x(m,n) 1S modeled, as in Chapter II, as 
a homogeneous random field with known mean and autocorre- 
lation function and it is considered independent of the 
target texture. The target region 1S Similarly modeled, 
but the amount of "a priori" knowledge, for most applica- 
tions, may be restricted to the mean and/or variance. 


The background region is modeled by 


zn) = Tr nl + v (m,n) (6.2) 
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Conditional probability density function (case 1) 
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The term ax(m,n) in (6.1) takes into account a possible 
correlation between the two regions. The weight a is a 
factor to be adjusted according to the specific nature of 
the picture. In the visible spectrum it seems reasonable 
to assume a replacement process, therefore a= 0. In the 
infrared spectrum it may assume a constant value at the 
edges of the target and zero inside this region, provided 


we assume that the sensor has rejected the background. 


me LIKELIHOOD RATIO TEST 

In this section we will construct a likelihood ratio to 
test the two hypotheses against a threshold. We will first 
consider the general case of deciding between two random 
fields, target and background, where the required "a priori" 
knowledge are the mean and autocorrelation function of the 
background and the mean and variance of the target and the 
measurement noise. After, we will work some special cases 
and derive its R.O.C. (Receiver Operating Characteristic). 
The detector developed here will be a Newman-Pearson detector 
[15], where the threshold is obtained from the R.O.C. sub- 
jected to the constraint that the probability of false 
alarm is less than some desired value. 

The two hypotheses are: 

Hypothesis H_ - The measurement y(m,n) is in the back- 


ground region: 
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Hypothesis Hy - The measurement y(m,n) is in the target 


region: 
VA = Ot v 


The decision rule will be applied to each pixel using 
the observation y(m,n) and the prediction of the pixel gray 


“a 


level Xo (m,n). Assume that the picture is being scanned 

row by row, from top to the bottom. For the sake of clarity, 
assume that the past pixels were correctly identified as 
belonging to the background region. The image filters 
developed before (recursive or hybrid) are used to give the 
prediction x (m,n). We are looking for a decision rule 


based only on y(m,n) and L 


Let's define the likelihood ratio: 


L(y,x,) = Io e (6.3) 
P (x,y | Ho? 
where p(x yi) is the joint probability density function, 


conditioned to the hypothesis Ho Or H)- 


Applying Bayes rule: 
p(x y | H) = p(x |Hj)P(y|x Hi) i = 0,1 


Since the prediction x, 1s independent of the hypothesis: 
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Hosting these results into (6.3): 


A DE 
Dyz) = In En. (6.4) 
p (y |x Eo? 


At this point we will make the assumption that the 
conditional probabilities density functions in (6.4) are 
Gaussian. 


We can write 


Ke 
Hi 
ll 
Y 
+ 
< 
| 
(D 


E 
HE 
pa 
ll 
Q 
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O 
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< 
i 
Q 
(D 
q 
3 


Since x, is known, the conditional densities are func- 
tions of (ve) and (v-ae tT), where So 1s the prediction 
ror. In figure (6.2) the histogram of se 1s shown for 
the noisy image in figure (5.8), using the first-step 
predictor of the two-dimensional recursive filter [3]. 
This result shows a curve quite close to a normal distri- 
bution, and, therefore, validates the above assumption for 
p(y |x H) and also for a provided the target 1s 
deterministic. 

Since the conditional probabilities density functions 
are assumed Gaussian, we need only to find its mean and 
Variance. Before proceeding, let's first define some 


IM 


Statistics for y and Ol 


I5I 








Figure 6.2 


Histogram of (ve) for the image in figure (5.8) 
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(a) Means 


To make the formulas more general, we consider here 


that the random fields x(m,n) and T(m,n) are not zero-mean. 


Since the fields are homogeneous, the mean is a constant 


and the same for all pixels. All the previous equations of 


the image filters remain valid, provided we change x to 
(x- x) and y to (y - x), where x is the mean value of the 
Field x(m,n). 


We will use the following means: 


(Gas) 


EI = E 

A T 

E(v) = 0 

E(y|H_) = Yo = x 
ERR =, = ox UT 
E (x) = E(x) = x 


(b) Variances 


The variances are given by 
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V(t) = Im 
Viiv = -f 
(6.6) 
Viy|H.) = o A 
O y 
O 
= 2 2 2 2 
V(y|H,) = y. Nor RC Wer 
(c) Covariances 
We have the following covariances 
Cov(v,X) = Cov (v, Xp) = 0 
Cov (v,y;) = Y, 1 = 0,1 
Cov (x, rx) = b (6027) 
Cov (xp Yo? = Cov (x, x) = b 
SNAP ASES ES = ab 
(d) Coefficients of Correlation 
The coefficients of correlation are 
A E Gentile e N bo 
p 50 m 
Cov (x_,y_) b 
Po = = (6:8) 
Geo g O 
Es P “e 


154 
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PY, P Y] 





Now, returning to equation (6.4) and using the Gaussian 


pdf's, we have 


A 2 “ 2 
A [y-E(y|x_H_)] [y-E(y|x_H.)] 
i een A DI. 76,9) 


viy |x Eo] Vly[x_H,] 
where we have dropped some constants, since they are not 


important in the test against a threshold. 


The regression of the mean of y on x_ is 


p 
- a fe ee | 
Ely|% E,1 = Yi + a p,(x =x) i= 0,1 (6.10) 


substituting equations (6.5) and (6.8) into (6.10) gives 


Ely|x_H = 1-y)x + : 
[y | E al (1-y)x + YX, 
(6:11) 
Ely|x H] = a (l-y)x+T+ayx 
where we have defined: 
= b BR 
Fa + > o, (6.12) 
pP 


Having the means of the conditional pdf's, we find 


its variances: 
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Vly |x HG] 


giving 
Viy |x Ho] 


Using equations 


A 


COV 
P 


y 
ly lx, 


Vly]x_H,] 


and 


V ly |x HI 


(6. 


rX) 


Ho] 


2 2 

E{ fy -E(y|x,H,) | } 
ur = 2 
E{ [x+v- (l-y) x-yx_] 7} 


P 


ee 08 +v]?] 


2 Pina 2: E 
g7 + o O AL CON AA 
Y p Y ( p 


7) and (6.12) we have 


2 
= vg 
í P 
= q? - YO é + r (Scl) 
P 
El [y -E (y|x ón 
l pol 
ee 
E{ [ax+T+v-a (1-y) x-T-ayx,] } 
D = A a 

E{ [a (x=x) + (T-T) -ay (x 7X) +v] Í 

222 po. 2 2 
A gos g y So + Sm + r (6.14) 
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Substituting equations (6.11), (6.13) and (6.14) into 


(6.9) gives 


— ^ 12 — — Se 
a E ] [y-a(1-y)x-T-ayx_] 
L(y,x_) = ES. 
p 2 Dit 2 2 
Ore OO ke DAI O t9- FE 
Y % Een 
Let's define: 
M = (l-y)x 
Homo.) 
2 
N = o ~yo 
Oy. p 
2 2 A — 2 
A (Y=Y 8 TM) (y-ayx,,-T-aM) 
L x = - 6.16 
(Yı p? N+r an eee 2 (6 
Ir 
2 


(a*Ntr+0q,°) (y¥-YX,-M) “= (Nex) (y-ayx -T-aM) 
= 2 


(Ntr) (a°Ntr+o, ) 
The denominator iS a positive constant, because 


2 22 2 
N = gog = O = o (l- 
Y D ( 2 


2) > 0 


Dropping the denominator and simplifying the numerator 


of (6.16), the likelihood ratio becomes 


Nyx) = [(y-H) -y (x,-#) 1° E S [ya RT) may (8) 1 


(Gre) 
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where 


Ze 2 
N+r lo. Chou Eur 


S a a (02.18) 
2 2 2 Zine 2 
a N+tr+om a one )o +r+o„ 


Pa 


Given the prediction re and the measurement y, the 
decision rule consists in computing the likelihood ratio 
L(x, ry) and compares it against a threhold Ny if greater, 


we choose hypothesis Hj + if not we choose hypothesis Ho: 


a ... Hypothesis Hy 


en ... Hypothesis AS 


En figure 6.3 a block diagram of the estimator-detector 
ilter is A It is assumed that the background pre- 
dominates in the picture and the objective is the detection 
of a smaller texture that we call target. The background 
texture is modeled as in chapter II. To accomplish such 
modeling, however, a picture without target is needed. In 
an actual application it is realistic to assume that 
we have a typical realization of this process (background) 
before the target enters in the scene. For 
example, assume that we are dealing with a surveillance 
problem, where some area has been scanned many times looking 


for a very low contrast target. Therefore a recursive or 


hybrid filter can be designed to predict and estimate the 
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background. About the target we need only to have a good 
guess about the mean and variance of its gray level. In 
cases of very low contrast the mean may be assumed equal 
to that of the background. The variance, in many cases, 
may be assumed much smaller than that of the background. 
Also, in many cases, the observation noise is much smaller 
than the background. 

The estimator-detector filter of figure 6.3 scans the 
image row by row from top to bottom. Assume that the 
previous pixels were correctly decided as belonging to 
the background. The decision box receives the measurement 
y and the prediction x and computes the likelihood ratio 
L(y, x) from equation (6.17). If this ratio is less than a 
chosen threshold n, the pixel (m,n) is considered background 
and the observation y is used to find the background estimate 
inn). If the ratio is greater than n, the pixel is con- 
Sidered as belonging to the target region. While in the 
target region the background filter cannot use the observa- 
tion y, unless we have "a priori" knowledge about the target 
gray level. Since we are considering that the target region 
is small, the background filter may run in its prediction 
mode without much degradation, provided the dynamic model 
ls accurate enough. Observe that both recursive and hybrid 
predictors can be used, the only difference being the box 
named "improved prediction" that applies only to the hybrid 


predictors. Since the three-dimensional predictors 


presented in chapter V, makes use of the measurements in 
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the previous frame, we can expect a good performance of 


those filters in the target region. 


The gray level T(m,n) of the target texture can be 
estimated using the background prediction Xo and the 


observation y|H,. 


AN 


T(m,n) = kY} + K 5X, 


where 


v= y|H, = OX + T tv 


The estimation error is 


e = T-T = T- k (ax+T+v) - KXp (6.19) 


k = tit, (ax+T) ] (6.20) 
X 


Using equations (6.19) and (6.20) the weights Ky and 


k, can be found using the minimum-mean-squared error criteria: 


aXT(o 2-b) + zo + Ro ©) 
C oS 5 > mn) = =a ee (6.21) 
a x (og Le SAB ja Aie e -b)+x In +T So +x r 
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Observe that if we neglect the observation noise 


(r = 0) and the prediction error (x_ = x, b = Ze = oh 


equations (6.21) and (6.20) become 


thus: 
T = Y = ox ED) 


Equation (6.22) is intuitively appealing, since it is 
exactly what we expect in the absence of observation noise 
and with an ideal predictor. The block diagram of figure 


(6.3) considers this simplified case. 


D. SPECIAL CASES 
In this section we will study some particular cases of 
detection. The detector developed in the previous section 
will be applied to simple cases of practical interest and 
its performance (R.O.C.) will be derived. 
2 


l. Case I: a=], Im = 0 


In this case the target gray level is deterministic 
and immersed in the background. The target gray level 
T is an unknown constant and we also don't know the target 
Shape or location in the picture. The objective is to 
Suppress the background in order to enhance the target. 


As before, the recursive or hubrid filters are used to 
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predict the background x and the detector, using this 
prediction and the observation y, decide whether the pixel 
(m,n) is target or background. The output is a binary 
image where l is target and 0 is background. 


The likelihood ratio, equation (6.17), becomes 


ra E ee? Ea 
L{y,x,) = [(y-x) (x, x) ] ot ly=x—-T) Y xo x) ] 


where 


Thus 
ne = 2T[(y-x)-y(x_-x)] - T 
(y Di [ (y )= Y ye) 
Dropping the constants, this ratio becomes: 


EX es) => (VaR = we (6-25) 


where 


un 
Y Ae 


D P 


Equation (6.23) is intuitively appealing, because it is 


exactly the residue of the recursive filters when y = 1. 
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This would happen if the filters were optimum. Such a 
decision rule with y = 1 was used in [5] as a generalization 
of the one-dimensional Kalman filtering technique. However, 
we have shown here that such an extension is not correct, 
because the filters are not optimum, and, consequently, 
the residue (or innovation process [16]) does not correspond 
to the likelihood ratio. However, it is a good approxima- 
tion to make y = 1 for the recursive filters, in cases of 
low observation noise and highly correlated background. It 
can be verified that y = 1 for the hybrid filters. 

Now, we will make the assumption that y and E have 
a joint Gaussian distribution in order to find the performance 
of the detector. Previously we had assumed that the con- 
ditional pdf's were normally distributed and we validated 
Such assumption with experimental results, but this does 
not mean that the joint distribution is also Gaussian. 
However, it is analytically convenient here, in order to 
have simple results. 

If y and x, are jointly Gaussian, the likelihood 


ratio L is also Gaussian and its pdf is 


— 2 
- (2-28) 
2, (2) E E exp (6.24) 
A Sr 2 oy 
where: 
MI i = 0,1 
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2 
Op = 0 + yO NOS (625) 


li 
m 
| 
O 


The variable L has the same variance for both hypothesis 


E and H 


A ¡+ Only the means are different as giyen by 


E(L|H,) 


i 
| 
o” 
i 
$ 
+ 
= 
| 
x | 
i 
© 


EN ES R (SETEV) er 


In figure 6.4 the conditional pdf's of the likeli- 
hood ratio L(Y: Xp) is shown. 


The probability of false alarm is 


FU 
ti 


prob (L|H, > n) 


CO 


L f 2 2 
= = exp(-2£ /2o0. )dg 
/27 de, n L 
É i exp (-u“/2) du 
Y 2T 
n/a, 
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PR = erfc(n/o,) (6.26) 


The probability of detection is 


da Ss prob (L|H, > n) 
= —t— Jj expr-(2-7)?/20,2142 
Y 27 Or 
n 
Pp. = erfe(l>) cz) 
g 
L 
where 
l 2 
erfc(a) = — f exp(-u“/2)du 
Y 21 g 


The choice of the threshold can be done, for 
example, by specifying the maximum admissible probability 
of false alarm; in this case we have aNeyman-Pearson detector 
[15]. Another criteria is to minimize the overall probability 


BE error: 


FU 
il 


HG A e P(H,)P(E|H,) 


P (Ho) PF + P(H,) (1-P) 
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If we know P(H_) and P(H,), an optimum threshold 
can be found that minimizes Poe If this is not the case, 
the best we can do is to assume that the hypotheses are 
equally likely, P(H_) = P(H,) =e. With this choice, 


it can be shown that the optimum threshold is n = T/2 and 


Pa 7 Pp = 1 ~- Poe 
The implementation of the estimator-detector for 
this case is similar to the one shown in figure 6.3. In 


this case, since we assume a target with constant gray 
level, and, also assuming neglectable observation noise, 

we might use the measurement in the recursive estimation, 
while the filter is in the target region, provided that T 
ES subtracted from y|H,- The best estimate of T must occur 


at the edges of the target region, because there, the 


background prediction Xo is the best. 


2. Case II: Q = 1l, r Á 0 


In this case the target gray level is a random 
variable with known mean T and variance a and, as before, 
we don't know its shape or location in the picture. 

The likelihood ratio, equation (6.17), becomes 


[ (y=) -Y (ex) ee (xx) 7]? 


E(y, X) 


2 


(1-6) (y 7 12 + 28T[(y-X) =y (4,5) 1-87 


(0.28) 
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where 


2-2 


zes Korn r 
AA 22 
2 2 2 
(1 > XO HE+, 
and 
oi 
1-6 = 23 o e L 
(1-9, Jo HEO m 


Rearranging equation (6.28) in order to have a 


quadratic form, and observing that (1-6) is positive: 


L(y,x,) A E 
nn [ (yx) Y (x 7) ] t 281 (y=x)=Y (xx) ] 


ce = eee 


[ (Y-F) -y (x 7x) +6] - 88 + T) 


Dropping the constants, this ratio becomes 


Liy,x.) = [y-E)-y(x -R)+8]? (6.29) 
p p 
where 

2 g? r 

B = DEZE ) u sa! 
T dh 
Y = e 
a D 
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To calculate the performance of the detector, we 
make again the assumption that y and Xo have a joint 
Gaussian distribution. 


Define the random variable Z: 
Z = (y-x) - y (X7X) + 8 


Since y and Xo have a bivariate normal distribution, 


Z is also normally distributed and its pdf is 


1 = (2-7)? 
Pz (Z) = — exp I —] (6:30) 
VETO 210) 
Z Z 
z = y - x + 8 
o é = g + y 0 = 20. Y a LO] 
Z Y Pp Y P 
2 Dis dd 2 22 
=p 0 = g = gq - o 
A E ete aD 
Since L = 2°, TES pdf 1s 
P(2) = + (P, (va) + P,(-/2)1 
2/4 
Thus 
1 CATA Se 
P, (2) = ooo! + o 
20 20 
207/22 Z Z (6.31) 
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The mean and variance of Los can be computed 


using the Moment Generating Function of 2. 


m Ze = — VA 
M, (t) = E(e“~) = exp[tz + 50, | 
E(L) = E(22) = MMO) = 2° +0,” 
V(L) = E(L?) - E?(L) = Eizd) - 22(22) 

E IV B Ii 2 

O [M(0)] 

Therefore 
E(L) = Zº 4 u 
(6.32) 

V(L) = 20.%(0.% + 22°) 


A 


The distance between the means of Byers) under 


the hypothesis Ho and H is calculated using (6.30) and 


A 
M0. 32): 

2 | = 
d = E(LjH,) E(L|H ) 
— In 2 lia E 

ECLIH = T+ +0 +0 +r+Yy 0_ -2p.0._ OU 
(L|H,) (T+8) p tFtY"O, 2019, 9, 

= (THB) “to eo rm yo À 





E(L|H) e + y%0 wo Sp 
cece o o é 
p 
Baur 
p 
thus: 
d = T + 29T + Sm (6.33) 


From equation (6.33) we can see that the distance 


between the means 1S non-zero even for a target with zero 


mean. 
The probability of detection is 
Py = prob(2|H, > n) = fp, (2/H,)4d2 . 
n 
Using equation (6.31): 
E -OT -= z)" 
ii ak il 
Dan f exp [ 5 142 
5, 2m 20 
à i 1 
co — ,2 
1 1 (vê + 24) 
en 5 ji — exp a a jde (6.34) 
Gian ei 20 
dl Aa 2 
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Using the complementary error function (erfc), 
equation (6.34) becomes 


E = erfel—(n/? -2,)1+ erfcl— (n + z7)] 


E oa 





(6255) 


Similarly, the probability of false alarm is given 


by 


CO 








= | = | 
Rs cont a a 
n 
(6.36) 
2 RE to a 
Pp = erfc[= (n z 11 + erfcl- (n tZ) 
Z Z 
O O 
The values of z and Or are given by 
= = r 2 0° 
T T 
a 2 2 
ce = (1-0 5 jo + In +r 
1 
2 
a ee 2 PR CAR 
u N Po o o 
T T 
O e iaa o + Y 
Z p 


To analyze the detector in terms of some meaningful 


parameters, let's define the ratios: 
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2 


BNR = —- (background-to-noise ratio) 
T +0, 
TBR = 3 —(target-to-background ratio) 
O 
=2 
Ep 
RTT = E 5 i 
T 


Normalizing the equations with respect to the 


Background (0° = 1), equations (6.35) and (6.36) become 
= n = BART 
Py = erfcl = TER b+RTT] 
1+RTT 


(Ga 39) 


de sro) em + Yb+RTT] 
AFTIRTT 


E ect - /b] + extol |? + yb] (6.40) 


HU 
tt 


where 


Bu muB 
a = att “p ) 
(6.41) 
a RTT 
Br a TBE (L + RTT) 


Using (6.39) and (6.40) we have calculated the 
R.O.C. (Receiver Operating Characteristic) of the detector. 


In figure 6.5 the effect of the target-to-background ratio 
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Figure 6.5 


Curve of Performance changing the 
target-to-background ratio (TBR) 
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(TBR) is shown. As expected, the detector performance is 
highly dependent of this ratio. For a probability of 

false alarm of 10% the probability of detection varies 

from 67% to 97% when TBR varies from -3 to 3 dB. In 

figure (6.6) the effect of the background-to-noise ratio 

is shown. Since this ratio is responsible for the accuracy 
of the background prediction E it is also an important 
factor in the detection. From this figure we can see that, 
for a probability of false alarm of 10%, the probability of 
detection increased from 63% to 93% when the BNR varied 
Brom 10 to 30 dB. 

In figure (6.7) the effect of the randomness of the 
target gray level, represented by the ratio RTT is shown. 
We can also observe a considerable improvement in the 
detection when this ratio is increased. 


3. Case III: a= 0, Im = 0 


In this case the target gray level is deterministic 
With known value, though its shape and location is unknown. 
The target and background regions are completely uncorrelated 
(a = 0), thus it is the case, for example, of images in 
the visible spectrum where there 1s target or background, 
but not the addition of these textures. 

Tmhe likelinocod ratio, equation (6.17), for this 


case, becomes 


AN ~ 


Liy,x.) = [ (y-X) -y (x 77) 1“ - öly-T]? (6.42) 


1575 





o 





Figure 6.6 


Curves of performance changing the 
background-to-noise ratio (BNR) 
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Figure 6.7 


Curves of performance changing the ratio RTT 
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Let's find the distance between the means, under 


the two hypotheses Ho and H}: 


A 


Dex) © (y-¥) +y? (x =X) "= 2y (x 7E) (y-F) -6 (y-1) 7 
VAE = x+v 
E(LIH) = 0° +rty°o Ÿ-2Yb-6E [ (y-x)-(T-X) 1 

= 4e vo 7-6 (0° +7) -6 (7-x) * 
E(LIH ) = (1-6) (07+r) - yo, - 6 (T-x)* (6.43) 
A EO RR 
E(L|H,) = Grea cee o ve,” (6.44) 


Subtracting (6.43) from (6.44) we find the distance 


a= ELEMENT | EH") 


(6.45) 


2 2 


Ou 
il 


(o Cr) É + (26, +8-1)0 


Observe that the distance is non-zero, even for tar- 
gets with the same gray level as the mean of the background 


E = x). 
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The performance of the detector can be found, 
provided we assume, as before, a joint Normal distribution 
for y and a however it is not possible to have an elegant 
analytical solution as in the other cases. 


4. Case IV: a= 0, © # 0 


This case might cover the model of some images with 
two quite distinct textures. A difficult problem, for 
example, is the restoration of pictures with two different 
regions, one is the sky and the other the sea. This is a 
typical case where we cannot model the image as a homogeneous 
random field. One approach to this problem is to model each 
texture, separately, as two homogeneous random fields, by 
the method of Chapter II, and design one image filter for 
each texture. To apply these filters, first we have to 
decide to which texture the pixel, to be processed, belongs. 
Such a decision can be made by the detector of the present 
case. In figure (6.8) the block diagram for this situation 
is shown. Following the same notation as before, we have 
a filter for the background x (one texture) and another 
Filter for the target T (the other texture). The threshold 
decision selects the proper filter, while the other filter 
runs in its prediction mode. 

The likelihood ratio for this case, equation (6.17), 


is 
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Let's compute the distance between the means: 


E(L|H_) = (1-8) (07 +r) -y"0 2-8 (T-X) 7 

air Eee ae (1-8) (op +r) = yo," + 52 

d = E(L|H,) - E(L|H ) 

d = (1#6)(T - x)% + (1-8) (0-0) + 2050? (6.47) 


Equation (6.47) shows clearly the differences 
between the two textures (mean and variance). 

As in Case III, the performance of this detector 
can be found using the Gaussian hypothesis, however it is 


not possible to have an elegant analytical solution. 


180 















FILTER FOR 
TEATURE x 





DECISION 


L > 
ni 





> 






ETERERIEOR 
TEXIURE T 






Figure 6.8 


Estimator-detector filter for two textures 
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VII. TARGET TRACKING 


A. INTRODUCTION 
In the previous chapters we have developed statistical 
image filters that exploit spatial and/or temporal redun- 
dancies existing in the pictures. One application of those 
filters is the restoration of images contaminated with 
observation noise. In this application, the objective is 
to remove tne observation noise, in order to recover the 
original image. Another application is to extract some 
desired information contained in the pictures, which we have 
called target. In this case, we want to suppress the rest 
of the information, which we have called background, in 
order to enhance the target. 
The target detection and tracking problem, by means 
of pictures sequenced in time, may be divided into three 
phases: 
(a) Extraction of the targets from the background, 
creating a binary picture (l-target, O-background). 
(b) Recognition of the target(s) that we are interested 
tracking. 
(c) Track the target(s) from frame to frame. 
This research addressed the first and third phases. 
The recognition problem, for our present purpose, will be 
considered as solved by a human operator, or, in the case 


of an automatic system, by some existing algorithms [18], 
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[19] of pattern recognition. The extraction problem was 
covered by the previous chapters. 

In this chapter we address the tracking of targets 
from frame to frame. The basic idea is to model the target 
movement within the pictures, and to construct a Kalman 
filter to track the target centroid, using the image filters 
as the measurement device. Of course, others relevant 
points of the target might be tracked, in a similar way, 
in order to find, for example, the target attitude. The 
image filters referred to here are the three-dimensional 
recursive and hybrid predictors, presented in Chapters IV 
and V, respectively, and the detector developed in Chapter 


ET. 


B. TARGET TRACKING FILTER 

In this section we will construct a Kalman’filter to 
track the target centroid. Assume that we receive pictures 
at some specified rate, say 30 pictures/sec. Since the 
Pictures have only two dimensions, the target-image dynamics 
is related only to the components of the actual movement 
that are parallel to the picture-plane (see figure 7.1). 
The effect of the movement perpendicular to the picture is 
to increase or decrease the size of the target-image. If 
the target has the same aspect at all angles (for example: 
a sphere), then it is possible to compute such component , 


uSing successive frames. 
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The target-image becomes a point-target at large dis- 
tances, but, in general, it is a mass-target with some shape 
which changes when it maneuvers. Such change can have very 
interesting applications in the target tracking problem, 
since it gives useful information of target maneuvering. 

In what follows, we model the movement of the target 
centroid in the plane of the picture and construct a Kalman 
filter to estimate its x and y coordinates, and also to 


predict its position in the next frame. 






TARGET- IMAGE 


TARGET 





PICTURE-PLANE 





Figure 7.1 


Target dynamics 
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1. Dynamic Model 


Since we don't know the intention of the target 
driver, we will consider two random accelerations in the x 
and y directions, therefore the position of the target- 


image is given by 


aii) Mee ce) Seer) + Su, (x) 7 
(7.1) 

o 1 2 

y(k+1) = y(k) + y(k)T + $0, (k)T 


Equation (7.1) assumes that the time separation T 
between frames 1S small enough such that we can consider 
that the acceleration and velocities are constants between 
frames. 


The velocities are given by 


x (k+1) x (k) + w. (K)T 


2) 


Y KER) 


y(k) + w, (k)T 


where we have used in this derivation the simplified 
notation x(k) for xtkT). 
To use state variable notation, define the state 


mector x and the random forcing input w: 
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ES 
I 
E 
Il 


In state space notation, equations (7.1) and (7.2) 


become 


l T 0 0 
0 ] 0 0 
x(ktl) = x (k) 

0 0 I nT 
0 0 0 lái 

7/2 0 

dk 0 

+ > w(k) (es) 
0 7.72 
0 dE 


The random forcing input w is considered zero-mean 
white noise uncorrelated with x, and has the covariance 


Matrix below 


Elw(k)w (kK)] = Q = (7.4) 
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2. The Observation Vector 
The centroid of the target-image is measured by 
the image filters previously developed. The pixels are 
examined one by one and classified as target or background. 


After this classification the target centroid is computed 


by 
N 
= l u 
A ) n, Tim, ‚n,) E S = ) T(m n) 
k=1 k=1 
N 
Em 
mern m RS) 
k=1 


where (m, ,n,) are the coordinates of the pixels classified 
as target. If we create a binary picture with Tim, ,n,) = 1, 


equations (7.5) become 


N 
x = i in 
m N k 
k=1 
(760) 
= 2 57m 
Im N k 
k=1 
Since the image filters are not perfect, some 
target-pixels are missed and some background-pixels are 
taken as targets. Therefore the measurement of the centroid 


has an error: 
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where Vy and V, are the errors in the x and y coordinates, 
respectively. 

Let's analyze the errors Vy and Vo: First deter- 
mine if they areuncorrelated. In figure 7.2 assume that 
the image filter missed the target-pixels in region A. 

The result is an error in the observation of the centroid. 
Both x and y coordinates will be in error. In other situa- 
tions only one coordinate is wrong. In general, we can not 
draw any inference about one error by the knowledge of the 


other, therefore we will assume that Vi and v, are uncorrelated. 





TARGET 





Figure 7.2 


Measurement error 
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We must determine if the errors are correlated with: 
the coordinates (x,y). Since the background and the target 
gray levels are considered homogeneous random fields, 
there is no reason to think that the image filter will be 
more likely to make errors for some positions of the 
target within the frame, or at some special frames. On 
the other hand, the nature of the filter itself could be 
a source of correlation. Both image filters (recursive 
and hybrid) use recursion in space, therefore the background 
prediction is poor near the top left corner and also near 
the first row and column of all frames, because of the 
boundary transients of the filters. In the first frame the 
hybrid predictor does not make use of the time correlation 
in this frame. In the recursive predictor the first few 
frames have transients in recursion in time. The conclusion 
is that there are some positions (x,y) in the time-frame 
that have poor probability of detection. However, we have 
seen that the recursive filters converge very fast, there- 
fore, we will neglect such correlation, since it occurs only 
in a very small region. 

One must test to see if the errors are zero-mean. 
The detection error could happen at any pixel, but the bias 
is most likely in the direction that the picture is scanned 
(top left corner to right bottom corner), because when the 
filter is in the target region the background predictors, 


in general, can not use the observation y(m,n), which 
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resultsin a degradation in the detection. However, we 
expect that such degradation can be neglected for practical 
applications where the target has small size and the back- 
ground correlation is high. There are other schemes for pro- 
cessing that can be done to overcome this problem. One is 
for example, to process the pictures top-to-bottom and 
bottom-to-top, alternatively, from frame to frame. 

After these considerations we conclude that it is 
reasonable to consider the errors Ya and V, zero-mean white 
noise and uncorrelated with x and y. Such conclusions 
simplify considerably the following derivation and the 
result will be a very simple Kalman filter. 


Define the measurement and error vectors as below 
An Vy 
y = y = Ta) 


In state space notation, equation (7.7) becomes 
2) E x(k) + v(k) (7.9) 


where v(k) is a zero-mean stationary random process, uncor- 


related with x(k) and with the covariance matrix below 
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7 I 
IAS ES IÓ 
0 
3. Kalman Filter 
The dynamic model is 
ote = kxk) +T wk) 


where 


D", 


IX 
li 
E 


T91 


(2210) 


EA 


(Cale) 


(ESA) 


(la) 





The measurement model is given by 


y(k) = HF x(k) + v(k) (eLo) 
where 
1 0 0 0 
a (7.16) 
6 0 il: 0 
X m vi 
VA Mo Des) 
Ym Ya 


The random forcing input w(k) and the measurement 
noise WG are considered white noise uncorrelated between 
themselves and also with the states. Its covariance matrices 


are given by 


Or R = (7.18) 


The Kalman filter equations are [17]: 


x(k) = x'(k) + G(k)[y(k) - Hx'(k)] (7.19) 


t92 





INES ES cne estimated state vector and x' (k) is the 


predicted state vector given by 
x' (k) = 0 x(k-1) 2.20) 


G(k) is the filter gain, computed as below [17]: 


60) = P"(x) 4 (ep (00 +r]"? (7.21) 
P(k) = [I-G(k)HIP'(k) (7.22) 
P'(k+1) = 0P(kK)0* + ror? 7 23) 


where P(k) is the covariance matrix of the state error 
and P'(k) is the predicted covariance matrix of the state 
error. 


The estimated state error vector is 
EEE sto ek) (7.24) 
The covariance matrix of the state error is 


= 


P(k) = Ele(k)e' (k)] (7.25) 


The predicted state error vector is 


O) (22.26) 
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The predicted covariance matrix of the state error is 


P'(k) = Ele'(k)e'! 


(k) ] (O 
The gain G(k) is computed by initialization of 
equation (7.21) with P' (1) to have G(1), then find P(1) 
by equation (7.22) and P' (2) by equation (7.23). Following 
this procedure the matrices G(k), P(k) and P'(k) are 
calculated. 
Since the x and y coordinates of the target centroid, 
as well as its measurement errors Vy and Vos are uncorre- 


lated, we can initialize equation (7.21) with the predicted 


covarlance of the state error below 


(oe eee 0 0 
(PO aan 0 0 
P'(1) = 225) 
: ee 
0 P34 (1) Pig fh) 





With such initialization the matrices G(k) and 


P(k) have the forms: 
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Working out the equations 


the results are 


g} (k) 


g (k) 0 
g, (K) 0 
0 g, (K) 
0 9 4 (K) 
P.a (k) 0 
P,, (K) 0 
0 P.. (k) 
0 Pag (X) 
AO) 
Pj, (k) 
en 
Pi, (K) 


RAF PK) 


1 Je 


Bel) 


RA A P-k) 


2 33 


P34 (K) 


Rowe >. (Ki) 


2 33 
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Enrougn 


(7.29) 
0 
0 
ON 
Pa (k) 
P ya (K) 
(7.23) 
(CS) 





I 
Py, (k+l) 


1 
Pl, (ktl) 
Ps, (K#1) 


P43 (k+l) 


ß 
P,(k+l) 


P44 (5+1) 


(ke) =  [l-9, (K)]P], (x) 


(k) = [1-g, (kK)1P;, (k) 


(k) = P3 (K) = gJ (K) P}, (kK) 
732.) 
33 (K) = [1-g, (k) ] P3, (k) 
34 (K) = [1-g, (k) ] P34 (k) 
44 (*) = Pig (Kk) -g4 (kK) Pa, (k) 
2 TA 
Py y (k) +2TP, (kK) +T P54 (k)+7T 91 
: . Le 
P12 (K)+TP5,(k)+5T Q- 
P.. (k) +770 
22 a 
072333 
2 1.4 
P33(k)+2TP,, (k)+T P44 (K) +7T Q, 
ee 
Pog (K)+TP}4(K)+5T Q- 
P (k) +270 
44 2 
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Initializing the filter equation with the initial 


Prediction x'(1), we have the estimated state vector 


x(k) = x'(k) + 9, (k) [x, (x)-x" (x)] 
x(k) = x'(k) + go (k) [x (k)-x'(k)] 
(7.34) 
y(k) = y'(k) + 9, (k) [y (k)-y'! (k)] 
vik) = y'(k) + 9, (k) ly, (k)-y! (k)] 
The predicted state vector is 
x' (kt+1) = os + T os 
x'(k+1) = a 
(7 335) 
Vater) = y (k) + T 7 (k) 
yi (k+l) = y(k) 


In figure 7.3 we show a flowchart for implementation 


of the filter. 
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Figure 7.3 


, Kalman filter flowchart 
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C. FILTER APPLICATION 
The tracking filter developed in this chapter is used, 
primarily, to provide the sensor servo with the predicted 
centroid coordinates for the next frame. Besides this 
application, the information of the predicted position of 
the target can be used: 
(a) To reduce the image processing, 
(b) To improve the performance of the target detector. 
The reduction in the image processing can be achieved 
by delimiting to a small area to be processed, centered 
at the predicted position of the target centroid. Instead 
of processing, for example, a 100x100 picture we might 
work only in a 20x20 region, depending on the target size. 
With such reduction we could afford to make more processing 
in this smaller area, in order to improve the performance 
of the target detector. 
Some examples of this extra processing are: 
(a) Scan the picture in two directions, one starting 
at the top left corner and the other at the bottom 
right corner. The results are two binary pictures 
(l-target, O-background) with greater probability 
of error for the pixels located in the transient 
region of the recursion in space (see figure 7.4). 
One simple way to combine the two estimates is to 
pick the half of each picture opposite to the starting 
point, because, in this way we remove most of the 


unfavorable region of detection. 
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(b) 


If more processing is allowable one might scan the 
picture four times starting at each corner and pick 
the quarter of each estimated picture opposite to the 
starting point. In this case the detection is 
completely free of the transient region. 

In the detection of the target edges the direction 

of scanning the picture is very important. Since a 
good prediction of the background is needed and, in 
general, such prediction is degraded while in the 
target region, the result is that the detection of the 
edge at the left is better than the one at right 

(see figure 7.5). One way to overcome this problem 
is to use the estimator-detector only in favorable 
conditions. Assume that the prediction of the target 
centroid is good. Process first the left hand part 
of the picture until the line AB (see figure 7.5). 
Next, process the right hand part, starting at the 
CEOP right corner. In this way, the filter will work 
in more favorable conditions, since it will face both 
left and right edges with good background predictions. 
The implementation of such a scheme may be simplified 
if we invert the columns of the right hand part before 
and after the procesSing. In this way, the filter 
estimates the left hand part and it is reinitialized 
at the boundary (line AB) to proceed to the estimation 
of the inverted right hand part. Therefore, the only 


additional process would be the required inversion. 
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Figure 7.5 


Edge Detection 
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VIII. COMPUTER SIMULATION 


em. DESCRIPTION OF THE SYSTEM 

A computer simulation of the target tracking and detection 
problem was implemented. The block diagram of the system is 
presented in figure 8.1. It is composed of three parts: 
time-frame generation, image processing and target tracking. 
The generator of the time-frame creates images at some speci- 
fied rate and characteristics. The images contain the target 


and are contaminated with additive white Gaussian noise. The 


target has a random movement from frame to frame. The images 
are processed in two modes: search and window. The search 
mode is used to acquire the target. In this mode the whole 


picture is processed in order to extract the target from the 
background. It is assumed that a recognition phase exists 
to select the desired target. In this simulation we have 
used only one target to avoid such a phase. After the 
acquisition phase, the images are processed only in a window 
centered at the predicted position of the target centroid. 
The output of the image filter is the measurement of the 
coordinates (4 ms Bm) of the target centroid. The target 
tracking filter receives these measurements and updates 


NA 


the estimation of the target centroid (a,ß) and also 
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computes the prediction for the next frame (01 Bo? which 
is feedback to the image filter. 

More details about each block of figure 8.1 are given 
in the rest of this section. 

1. Time-frame Generation 


The background has the autocorrelation function: 
R(i,j,K) = 0 O 


The ¿dynamic model that generates a random field 


with such an autocorrelation function is 


x(m,n,t) = A x + w(m,n,t) 
where: 
A = [P Phet “PLP TPL TOO PPP] 
x" = f[x(m-l,n,t) x(m,n-l,t) x(m,n,t-1) x(m-1,n-1,t) 
x(m-1,n,t-1) x(m,n-1,t-1) x(m-1,n-1,t-1)] 
2 _ _ 2 _ 2 y 2 u 2 
E(w (m,n,t)] = Q = a (1 O, er Ph CE Pr ) 


The block "Background generator" implements this 
dynamic model by driving the equation with zero-mean white 


Gaussian noise and variance Q. 
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The block “Target Generator" creates a binary image 
(T - target, 0 - background). The target has an arbitrary 
shape and it is centered at the location (x,y) given by the 
block "centroid movement". Such movement is generated using 
the dynamic model of equation (7.1). 

The block "noise generator" creates an image with 
zero-mean white GausSian noise at some specified variance. 

Then, those three images are added in order to 


Simulate the situation addressed in case I of chapter VI. 


SS SY 


y | Hz x + T+y 


where: 
x - background gray level 
v = observation noise 
T - target gray level 


2. Image Processing 


The three-dimensional recursive and hybrid filters 
were implemented in order to predict the background gray 


level x_. The decision rule using x, and the observation 


y is 
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A > n -... Hypothesis Hy 


<n .... Hypothesis Ay 


The block diagram of the combined estimator-detector 
filter is presented in figure 6.3, particularized to case 
mi(a= l, oa aoe One. 


The constant y was defined (see equations 6.7 


End 6.12) by 


y `= T 
= 2 E 2 as p 
P p 
where: 
2 aa A 
g = E = E + 
e on Lx ep)! 


COV X,X ro x e 
( E! es! 


For the recursive filter the prediction error en 
is not orthogonal to Kr because the filter is not optimum. 
However, the hybrid filter is such that the prediction error 


eo is orthogonal to y therefore: 
l ee for the hybrid filter 


for the recursive filter 
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The threshold n is chosen such that the probability 
of a false alarm is below some desired value. From equation 


(6.26) 


In the simulation the threshold can be specified by 


the operator or automatically set at the previous value 


E = 30, > ¡SES + r (see equation 6.25) 


Observe that this criteria is completely independent 
of "a priori" knowledge of the target gray level. 


We have defined the following ratios: 





TBR, = Target-to-background ratio at input 
TBR_ = Target-to-background ratio at output 
BNR = Background-to-noise ratio 
TBR, 
EG. = TBR, Processing gain 
2 
TBR. = = CD 
als 2 
DEE 
2 
E” [L|H,] a 
TBR = —— M = — 8.2) 
O 2 2 3 2 
L L 
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o +r BNR + 1 


E (1 - p,°) BNR +1 





BNR = 


From (8.3) we can see that one upper bound for 


P.G. is 


P.G. < BNR + 1 


Therefore, the observation noise has to be well 
below the background in order to have a substantial target 
enhancement. However, in many applications, the observation 
noise is really neglectable when compared with the back- 
ground and, therefore, substantial target enhancement (or 
background suppression) can be achieved with such a method. 

After the detection of each target-pixel, the 


centroid is computed by 


where (m,n, ) are the coordinates of the 1 target-pixel. 
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Cee alee Tracking Filter 


The target tracking filter is the implementation of 
the Kalman filter constructed in Chapter VII. It receives 


the target centroid measurement (%_,8_) from the image filter 


m’ m 

and computes the estimated position (0,8), as well as the 
predicted position in the next frame o The filter 
gains are computed on-line according to the flowchart of 


Elgure 7.3. 


Er SELECTED RESULTS 

In this section we present some relevant results of 
the simulation. The pictures have the size 34x70, compati- 
ble with the Tektronics 4012 Display. The images are 
pictorially displayed with 8 levels of quantization obtained 
by superimposing some characters. The target has a 
"diamond" shape, but it can be easily changed to other 
shapes enclosed by a 7x7 matrix. In the tracking mode 
the image processing is only applied to a 20x20 window 
centered at the predicted centroid of the target. 

l. Background Prediction 

Since the decision rule is highly dependent on the 

background prediction, Jenne measured the prediction 
error for several situations and compared with the theoreti- 
cal results. The results are presented in tables 8.1 
through 8.3. Observe that the theoretical and experimental 
results have the same order of magnitude. We can also 


observe that the experimental results are always greater 
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TABLE 8.1 





SDSHYERTE 





3D-RECURSIVE 











NOISE PREDICTOR PREDICTOR 
VARIANCE 
THEORY EXPERIMENT THEORY EXPERIMENT 
"293 «306 282 -360 


-2958 "332 






PICTURE VARIANCE = 1.0 


TABLE 8.2 


SP=RECURSIVE Span Ven. ob 
NOISE PREDICTOR PREDICTOR 


VARIANCE 
THEORY EXPERIMENT THEORY EXPERIMENT 
E 291 . 364 -305 >55 
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TABLE 8.3 







3D-RECURSIVE SD=-HYBERID 














NOISE PREDICTOR PREDICTOR 

VARIANCE 
338 3478 
217 52 
292 >30 
280 295 
.216 „244 

than the theoretical ones. Such discrepancy is due to the 


fact that the theoretical values are computed for the 
Steady-state condition, but the experimental values are 

the result of averaging over a small ensemble (20 frames) 

and small images (34x70), which means a reasonable influence 
of the transient of the image filters. 

Therefore, we consider that these results are a 
reasonable validation for the theoretical methods of filter 
design and performance evaluation presented before. 

2. Transition Region 
To observe the transient of the image filters, we 


have applied the decision rule to a time sequence of images 
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(time-frame) without target. The threshold was set at 
n = 0.35. The characteristics of the images are: 
2 


o = 1, O.) = OD = 0.93, 


E = 0.95, BNR = 30 dB 


Pt 
The results for the recursive filter are shown in figures 
8.2 through 8.5. As mentioned before, the 3D-recursive 
filter gives degraded predictions at the first few frames 
and near the first row and column for all frames, due to 

the time and space recursion, reSpectively. Observe that 

in frame 2 (figure 8.2) the false alarms occurred mostly 
near the first column and row, but they can also be Seen 

in the rest of the picture. In frame 10 (figure 8.4) the 
Situation is better, becauSe it the Steady-state was reached 
in time, but the false alarms persist due to the transient 
in space. In chapter VII we have suggested some methods 

to reduce such undesirable effects of the transition region. 
In this simulation we have used a simple method that is 


to increase the threshold for the pixels at the first 5 


rows and columns. The results of such a method are shown 
in figures 8.3 and 8.5. Observe the substantial improvement 
in both frames 2 and 10. The false alarms of frame 2 


(figure 8.3) are basically due to the transient in time. 
The percentage of false alarm after frame 5, as in frame 
10 (figure 8.5), is quite close to the theoretical value. 
Of course, the probability of detection is reduced in the 
transient region, due to the higher threhold, however this 


region is very small for the usual size of pictures (100x100). 
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Figure 8.2 False alarms of the recursive 
filter at frame 2, using 
one threshold. 
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Figure 8.3 False alarms of the recursive 
filter at frame 2, using 
two thresholds. 
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Figure 8.4 False alarms of the recursive 
Filter at frame 10, using 
one threshold. 
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Figure 8.5 False alarms of the recursive 
filter at frame 10, using 
two thresholds. 
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The results for the hybrid filter are shown in 
figures 8.6 through 8.9. Since this filter does not use 
recursion in time, its transient is only in space. Simi- 
lar improvement was achieved by using two threholds, as 
can be seen in figures 8.7 and 8.9, compared against figures 
8.6 and 8.8. Observe that frames 2 and 10 (figures 8.7 
and 8.9) have almost the same number of false alarms, as 
Should be expected, since there 1S no transient in time. 

Comparing these results we see that the recursive 
filter presented better performance than the hybrid filter. 
This result is according to the theoretical and experimental 
values of tables 8.1 through 8.3, where the recursive filter 
has better performance for high background-to-noise ratios, 
as is the case (BNR = 30 dB). 

3. Background Suppression 

In what follows we present some examples of target 
enhancement (or background suppression). Several parameters 
can be varied to evaluate the estimator-detectors performance. 
In these examples we have used the same kind of background, 
But different values of target-to-background ratios. 


The characteristics of the background are 


P = 0:95; BNR = 30 dB 


The correlations in space (041 Ph! have the same order of 
magnitude as those observed in chapter II for several real- 


life pictures. The background-to-noise ratio (BNR) was 
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Figure 8.6 False alarms of the hybrid 


filter at frame 2, using 
one threshold. 
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Figure 8.7 False alarms of the hybrid 


filter at frame 2, using 
two thresholds. 
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Figure 8.8 False alarms of the hybrid 
filter at frame 10, using 


one threshold. 
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chosen so high, because we are interested here in those 
situations where the observation noise is neglectable when 
compared with the background, which is the undesirable 
texture to be suppressed, in order o enhance the target. 
The performance of the filters will be given by 


the ratios defined in equations (8.1) through (8.3). 


Example 1 


Filter: Recursive 


Threshold = 0.5 
miguces: — SalsOn through 8.15 
In table 8.4 the performance of the filter at 
several frames is shown. Observe that at frame 1 the 3D- 
recursive filter reduces to a 2D-recursive filter which does 
not use the correlation in time. Comparing figure 8.11 


with 8.13 and 8.15 we can see the great improvement resulting 


from the exploitation of the correlation in time. 


Example 2 


Breer: Hybrid 


Threshold = 0.5 


Figures: §.16 and 8.17 
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Figure 8.10 Original image at frame 1 (TBR. = 0 dB) 
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Figure 8.11 Frame 1 after the recursive filter 
SS = 12 dB) 
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Figure 8.12 Original image at frame 6 (TBR, = 0 dB) 
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Figure 8.13 Frame 6 after the recursive filter 
(TBR, =o) 7.3 dB) 
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Figure 8.14 Original image at frame 10 (TBR, = 0 dB) 
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Figure 8.15 Frame 10 after the recursive filter 
(TBR, = 17.3 dB) 
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Figure 8.16 Frame 6 after the hybrid filter 


(TBR, = 15 dB) 
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Figure 8.17 Frame 10 after the hybrid filter 
(TBR, = 15 dB) 
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TABLE 8.1 


SD-RECURGLI VE, FILTER 










TBR. = 0 dB TBR; = -3 dB 
PSE example 1) E 3) 
2 ws 33 20 SU 
3 3 3 12736 16.6 
4 4 4 14.8 17-3 
5 .0 0 1803 1023 
7 “2 ae 15738 ses 
9 | N 26.3 19 
T3 7 e, TES 19.6 
dr 59 .9 16.4 19.4 
19 “2 2 r6-0 12920 
20 o o 17,3 20.23 | 
AVERAGE OVER 
FRAMES 6 av N 16.3 197 


THROUGH 20 


This example is similar to the previous one, but the 
hybrid instead of the recursive filter is used. The frame 
1 is not shown, since both filters reduce to a 2D-recursive 
filter for this frame. Comparing the results we can see 
that the recursive filter (example 1) presented better 
performance than the hybrid one. As expected, the eg 
filter had almost constant performance (15 dB of processing 


gain) for all frames, since it has no transient in time. 
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Example 3 


Filter: Recursive 
TBR, = -3 dB 
Threshold = 0.35 


Figures: 82 P32ehraugnz23,23 


The filter performance for other frames is shown 
in table 8.4, where the effect of the transition in time 


at the first 5 frames is clearly shown. 
Example 4 


Biveer: Hybrid 
TBR. = -3 dB 

1 
Threshold = 0.4 


Figures: 8.24 and 8.25 


This example is similar to the previous one, but 
using the hybrid instead of the recursive filter. The 
processing gain was around 12.5 dB for all frames. From 
figures 8.24 and 8.25 we can see that TBR; = -3 dB is about 
mie limit of detection for the hybrid filter. 

4. Target Tracking 

In what follows we present some examples of target 
detection and tracking from frame to frame. The background 
has the same characteristics as before. The first five 


frames are entirely processed in order to acquire the target 
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Figure 8.18 Original image at frame l (TBR, = -3 dB) 
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Figure 8.19 Frame 1 after the recursive filter 


(TBR, = 9 dB) 
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Figure 8.20 Original image at frame 6 (TBR; = -3 dB) 
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Figure 8.21 Frame 6 after the recursive filter 


(TBR, = 15.8 dB) 
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Figure 8.22 Original image at frame 10 
(TBR, = -3 dB) 
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Figure 8.23 Frame 10 after the recursive filter 
(TBR, LO Oman) 
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Figure 8.24 Frame 6 after the hybrid filter 


(TBR, = 12.5 dB) 
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Figure 8.25 Frame 10 after the hybrid filter 


(TBR, = 12.5 dB) 
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(search mode). The window mode is implemented at frame 6 
and, thereafter, the pictures are processed only in a 
small 20 x20 window centered at the target predicted 
position. In the examples we have created binary pictures 


to emphasize the performance of the filters. 


Example 5 


Filter: Recursive 
TBR; = 0 dB 
Threshold = 0.5 


Higures: 8.26 through 8.32 


Example 6 

Filter: Hybrid 

TBR. = 0 dB 

T 
Threshold = 0.5 
Figures 8.33 through 8.36 
Observe that the performance of the hybrid filter 

is inferior to that of the recursive filter in example 5. 
It can also be seen, in figure 8.36, that the tracking is 
better after the first 5 frames. The reason is that in the 


window mode the measurement error of the centroid position 


1s limited by the size of the window. 
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Figure 8.26 Original image at frame 6 
(TBR, = 0 dB) 
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Figure 8.27 Frame 6 after the recursive filter 
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Figure 8.28 Original image at frame 12 
(TBR, = 0 dB) 
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Figure 8.29 Frame 12 after the recursive filter 
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Figure 8.30 Original image at frame 18 
(TER, = 0 dB) 
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Figure 8.31 Frame 18 after the recursive filter 
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Figure 8.32 Target tracking with the 
recursive filter (example 5) 
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Figure 8.33 Frame 6 after the hybrid filter 
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Figure 8.34 Frame 12 after the hybrıd filter 
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Figure 8.35 Frame 18 after the hybrid filter 
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Figure 8.36 Target tracking with the 
hybrid filter (example 6) 


235 








Example 7 


Filter: Recursive 
TBR. = -3 dB 

JL 
Threshold = 0.4 


Figures: 8.37 through 8.43 


Example 8 
Filter: Hybrid 
TBR, = -3 dB 
Threshold = 0.4 
Figures: 8.44 through 8.47 
Comparing examples 7 and 8 we see that both filters 


were able to track the target, but the recursive filter is 


clearly superior to the hybrid one. 
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Figure 8.37 Original image at frame 6 (TBR. = -3 dB) 
1 


CLCLNL12212112JLLL2999332112923292923133123293823099011331231303033011JL ILLIA 
LLEENRLEIIIITFTISCELIIIIILLELSIITIIICIILLIDIIEIITIBIELLELLELIISSILCEII3I3EITLI 
NSSSN sSELILLLLNSCCLILLLSSSIIIIIIIIILLLIIDILZIII0III3LLLLLIIIICLII 3303323 
2, ELLLLLL“SSLULELUNNSS3ICLC3333322LLI00ILI00ODIIIILLNLLI9330.3300000000 
SSSSSSLIILLILLLLLLILUN NSS IICLLIILLLL3300LL303033I223L[.0[IC000032000000000 
LLEL\LIIIIIIILLLILIILLLLLIOSI33ILLLL3303ILLL3303332330C0328900390888800000 
LLLLNLINNIIIILLLLLIILLLNLIILILLL NSS LII2ILLL320232222230099988998388 00070 
SJILLLIDNNDANNIIIULADILLLLIIIILLU NSS SLLIILLLILSICLIIICI0I98800000888U073 20 
JALLLINDIOIDIIIIMINGILLLLILLILLLN \CNLLIILLLLLILLLLLLIOI0MII0080300 330 
JWILLIBBDNANSRDOADIDAIILLLILIIIIINLLLLIINIITNILILLIILI202N022200830800090 
JILNSIIIARAADTINDIDNIILLNSLLELLLILNCLLLIIIILIICIN N LLLLI000032399890800000 
JINSLIIMNNOSTDOROIIIILLNSLLLULLUILLLLLLIIILLILL : \CNSSCII3D0OILLL0000033320 
LL\»LI32I232208030I0IIIILN SS SSSSSLLNSCLNSLLLNLLNS: \LNNSL333232200000033023 


JIULLI2330900D0DIII09LIIU.NN encantar ecr ters: KAINIIIIDOBSDODDDIL I 
0093390008388000800399)0 Uns en cer rr er sr rc re cv LODD]1110080840011] 
209L.D0n9A08980000DII9InssNEnnanN ore cr cresc cr co LITITIIIIDIDIDIILLA 
ANIIID8SDE3AB3DDDDDJJILLLLLLLUENL. :ı 0::: 00 rer 0 He NILLLLLLAILIZIIILLL 
gaB009994ASASADDADDI)LLLLLCLCll : :: vv: : 0,00: 0 Hr: 0 ee LILIILLLLLILII3 ILL 
41948989391 499890992B91I01ITL 1370077 :::::: 00.0: 0 0 0 00 0 8 LILAIILLLLILIIILLL 
3nn08g94agn 99998088 909 209939393972 +: 0:00:01: 11: LOIODII2I0I0DNNNNZ2I 
4314192202999 1989990001197) :: cc: e LOULIILL 5800001 
S199048834448888888029222121233022J : "Bt" NJILLIILLIJ000001J1L 
338a9889490000088893021221221231J ce MEM: +: : \LLLILILIJIB8000J1L 
33938088888809908988119111122912139 : :'':: sss: - +: LENCIILLIDSBSSD IIL 
oooneaagaa9nDoennesenoooo0300000089 :::: | BEBBBER: :  L]1LL21L110988830]J1L 
203938388388000000820030333333231209: :::::::  @EBBB: :: : : : : 302I202323220983380 I IL 
83008003398000033303303333322020: :::: ums: :53003303000998988033L 
NDDODODODSBODNNDNNNDSDDIDIIII00DDDD: ::::::: Br 1: LIIIIIII0D0DBRSBBIIL 
0n23300833000003900I3332333L.2330008 :: :::::90:.0:20::::3932099309398849888333 2 
0090900088930000999099ILICCLrraaaadO css cr rr q0990090004888020L 
JVIDAAAASOOIAALIABBIALLLLLLAAADAGD ss secs ss [99900In0gaaaaganIIl 
32303339333900IIII0IILILLLLLLL212I0: ::: 00:0, 0:00 000001 IDDODBODSBBRBBAn 232 
3393330D0DRI0OGIIIIRDILIILLLLLL009: : : :::.: 90990000008840nILiL 
2293130099099239801122121211112900:' "ts: 333333I300908989003JL LL 


Figure 8.38 Frame 6 after the recursive filter 
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Figure 8.39 Original image at frame 12 (TBR, = -3 dB) 
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Figure 8.40 Frame 12 after the recursive filter 
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Figure 8.41 Original image at frame 18 (TBR, = -3 dB) 
TAJ 35135088030 0000099I HMC aca 
ii cest 2 * 2222120 0ILLLL2001LL311L11000 08 SRB 
LLN 2220.02 3333093310 0091 2311000000839 
Mar EOAIITOIIIIL LOL LIL IL II OIEA 
GIL css io site LILLO LITIO 
LIL III LLL II) IOMA ADO 
LUN) o LILLE ELL IOLLLILI NI RARAS 
LOL MM IDIJILLILLLLILIDLLTLLL. LITIO RAROS 
LIL 2209. : BERRA: :: |. JOJILLLLLLSLLLIJLLLLSL LI 1000naeessssEBascena 
LIL.) ERSA © OLLE O OLLA 1700 00898908988880€800 
LLL. >> AIBR JILLS SSLLSLLLLIILLLL> >` NELLLI220000088009010 11 
BEN ARE LANNANALENLLIDIDILLLELALILII3I13900090080993IIL I 
LIL o a E e eL LLNS SSSSSCLIIILLLLENNLIII2INIDDDDBECDOIDIL I 
Jd cc eS e LLL LLL 93333030000030000 323.1 
IL: rr LEN EN: :NELLIIIELLLLLOI00081900030000D00ILI 
ERR Pere re SENN Nt NLULLLILIIILINODIIONLIIIILIIIIIIIILe 
Stier out ge ge No beeen: NSNULULILILENINNNN00303399L33393391 3 
133 o LEN NNE NUUNLLIDIDIILIGDDIIDLIIILNCULLIIDAILLL 
733 E NN ESS SNESSLIOOIIILIDBDIIOLODOINILLLIII ILL 
Ju? Soe LEN  NLSSSLLNLIIDDOOILIDDDIIOIDDAILIIIIIIILLL 
149 2 LE NE LLLL3080333330003103000310 3333033333 


332323LI13C.323230D8DDDIIILLLLLENN  :SSS:CNSSSCIHIDDOIIILIIDITLLITLIIGILALIIIIILLI 
211112 LIL)1)00)0000L2JJJ1IJIJILNSNLLNS UNSDIIDDOIIILLIIDILILLIICNDLIIIILLLL 
LLLLLL LULLIII2I00IILIIIIIICLNLLIIUENIDLLLI08800IL.I2II0I2IUL.LIICLI200003333 
LU LLLLLIILLIIIIDIILIIILLIUNN SLULITLLILLI0088080CII00220 2320277200650 3330 
»LELLLLLIIIIIIIITIILLLIIIIUN N SCUILLILLI0000000I2I000330223I0 22270 222212 
ULLI II LIILLLELIILE ALU LCIIILLCLCII33333130390000030003083913I033113990D3L LL 
2333303103333333301313100393309333000333312303000003008809085000091339993311: 
33I2930231390033313033L11308333030093333013030808800088008 0000330393911 + 
L1239004080008011L3LL LC CLLLOTJTO0080D0111311102880000D888884300092]111111L;-: 
3309388888009 003333939tcCc333000008893323123030009323229033033039023233333L* . 
1137388888 0003331 3. 1333008888880 333823030000030008008 0000393190230T1> +» 
¿3001483388880 07I831 1930080888888 030332303800033300033003003993933393931> > + 
1330388883988 000333302390 0o0888882a8880003307800003930032309I9ILL 1. 333IL 


Figure 8.42 Frame 18 after the recursive filter 


239 





40 


0 - true position 
A - measured position 
x - estimated position 


30 
29 
10 
4 
9 
a 19 29 39 40 50 69 T9 


Figure 8.43 Target tracking with the 
recursive filter (example 7) 
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Figure 8.44 Frame 6 after the hybrid filter 
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Figure 8.45 Frame 12 after the hybrid filter 
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Figure 8.46 Frame 18 after the hybrid filter 
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0 - true position 
A - measured position 
x - estimated position 





Figure 8.47 Target tracking with the 
hybrid filter (example 8) 
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IX. CONCLUSIONS 


A. SUMMARY OF RESULTS 

The experiments with several real life pictures allow 
us to draw two important conclusions concerning image 
modeling. First, the hypothesis of an autocorrelation 
function with separable kernels in the vertical and hori- 
zontal directions worked quite well. This conclusion is 
very important, because the derivation of the dynamic model 
is substantially simplified by avoiding the difficult 
problem of two-dimensional spectral factorization. Second, 
the first-order model, used by most researchers, iS a poor 
approximation for pictures with few details, although it is 
a very good model for pictures with many details. The auto- 
correlation function, suggested in this research, seems to 
be a good choice, because it includes the first and second 
order models, as particular cases, and also permits a simple 
method of parameter identification which can be easily 
implemented. 

The performance evaluation of the sub-optimum two- 
dimensional recursive filters of Habibi's type [3-5] allow 
us to draw some interesting conclusions. Although sub- 
optimum, these filters are not so far from optimality. In 
the worst case, the error variance was around 15% greater 
than the one obtained by the optimum interpolator. 


Rosenfeld's filter [3] is the best and presented degradation 
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around 8%. Shachar's filter [5] has the simplest algorithm 
for gain calculation which uses a very good approximation, 
Since its performance was quite close to that of a similar 
filter suggested here which calculates the gains without 
approximation. 

The reason why Rosenfeld's filter [3] is the best is 
because the estimation error is orthogonal to the obser- 
vation of the pixel under estimation, although this does 
not happen for the rest of the data set. This conclusion 
leads directly to the hybrid filters introduced in this 
research. From the results shown in Table 5.1, it can be 
seen that the non-recursive filter [6], suggested by Bar- 
Yehoshua [6], although using only 9 observations, presents 
better performance than the optimum (in the sense of 
recursive filtering) interpolator, for correlations as hich 
as 0.94. The conclusion is that most of the information 
about the pixel gray level under estimation resides on its 
closest neighbors. For this reason, the hybrid filter 
presents the best performance, because it makes use of all 
the observations in the neighborhood and also those farther 
observations used by the recursive filter. Experiments 
with real life pictures (figures 5.12 and 5.13) show the 
ability of the hybrid filter in smoothing out observation 
noise in pictures with signal-to-noise ratio as low as 
-3 dB. The processing gain for this case was 11.1 dB, which 


compares favorably with previous results. 
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For images sequenced in time, the three-dimensional 
recurSive and hybrid filters were developed in order to 
exploit the correlation in time. The recursive filter is 
adequate for the cases where the scene is the same during 
several frames, otherwise the exploitation of the time 
correlation is not advantageous, due to the filter transient 
at the first few frames. Its structure is intuitively 
appealing, Since it reduces to the one-dimensional Kalman 
filter for the case of no spatial correlation, and to a two- 
dimensional recursive filter for the case of no correlation 
in time. Numerical results show that the use of time corre- 
lation can be quite advantageous for the case of time 
correlation greater than spatial correlation (error variance 
is reduced by 33%), but it is not so advantageous for the 
opposite case (error variance is reduced by 6%). The hybrid 
Filter is adequate for the case where the scene changes at 
every two frames, Since it makes use of only the previous 
and the present frames. Of course, it can also be designed 
for the case of a fixed scenario (see Chapter V) in which 
case it improves the performance of the recursive filter. 

To detect targets from cluttered background images a 
likelihood ratio was constructed which uses the observation 
of the pixel and the background prediction for the pixel 
gray level. A quite general situation was considered. 

The target and background are considered as two kinds of 


textures. The image is modeled as a weighted combination of 
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three additive components: target, background and obser- 
vation noise. Therefore, it can accommodate the case of 
infrared pictures and also visible pictures, where there is 

a replacement process, target or background. The prediction 
of the background is given by the recursive and hybrid filter 
or, also, by a non-recursive filter. The likelihood ratio 
was particularized for four special cases of interest, which 
had their performance analyzed. The likelihood ratio enhanced 
the target-to-background ratio. The target is extracted by 
threshold detection. The threshold is chosen in order to 
maintain the probability of false alarm below some desired 
value. The performance of the detection process is highly 
related to the background-to-observation noise ratio, since 
this ratio is responsible for the background prediction. 

The target detection and tracking problem was simulated 
uSing comptuer generated images with characteristics Similar 
to those of real life images used in the other experiments. 
Both three~dimensional recursive and hybrid predictors were 
used and compared. The likelihood ratio, using both filters, 
was able to detect a target at target-to-background ratios 
as low as -3 dB. The processing gain was around 19 dB, for 
the recursive filter, and around 12.5 dB for the hybrid 
filter. The recursive filter presented better performance, 
as should be expected from the results of Tables 5.3 through 
5.5, where its prediction is better than that of the hybrid 


filter, for the case of high background-to-noise ratio. 
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B. SUGGESTIONS FOR FUTURE RESEARCH 

Many pictures are clearly non-homogeneous, therefore 
the development of space-varying models is the natural 
way to proceed in order to improve the model suggested in 
this research. The recursive techniques are directly 
amenable to the analysis of space-varying models. The 
reliance on "a priori" information is a very important 
area. The method of parameter identification suggested 
here is applicable for the case of constant parameters for 
the whole frame. In the non-homogeneous case such identi- 
fication has to be accomplished pixel by pixel, in order 
to improve the robustness of this statistical approach 
in face of modeling errors. 

Due to the enormous computational load of the optimum 
multi-dimensional recursive filters, and also because the - 
observations far away carry negligible information about 
the pixel to be estimated, the natural way is to look for 
Sub-optimal recursive filters that require less computation. 

The hybrid filters suggested here can be equally applied 
with. other recursive filters, as, for example, the reduced 
dimension filters proposed by Panda and Kak [9] and Woods 
and Radiwan [11]. 

The decision rule suggested here is directed to the 
detection of targets pixel by pixel, in order to be indepen- 
dent ofthe target shape. The choice of the threshold was 


based on the probability of a false alarm of each isolated 
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f 


decision. Since a wrong decision necessarily affects the 
next, a question arises concerning how to vary the threshold 
in order to keep the false alarm probability below the 


desired value. 
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